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Image  Processing  for  Optical  Engineering  Applications 

Matthew  B.  Weppner  and  Matt  Young 

Electromagnetic  Technology  Division 
National  Bureau  of  Standards 
Boulder,  Colorado  80303 

This  Internal  Report  describes  a palette  of  computer  image  processing 
algorithms  for  picture  processing  and  optical  fiber  characteri zati on. 

Keywords:  computer  processing;  Fourier  optics;  g-profile;  Gaussian  beam;  image 

processing;  images;  optical  fiber  characteri zati on ; optical  processing;  single- 
mode fiber 

1.  Introduction 

1.1  Statement  of  Purpose  and  Organization 

This  Internal  Report  describes  the  development  and  testing  of  image  pro- 
cessing software  designed  for  optical  engineering  applications.  Image  pro- 
cessing functions  in  this  software  include  two-dimensional  Fourier  transforms, 
convolution,  noise  reduction,  multiple  image  resolutions,  and  low-level  image 
processing  functions.  The  software  also  contains  image  information  display 
tools  including  Gaussian  beam  and  g-profile  characteri zati on  for  optical  fiber 
measurements.  The  necessary  image  file  input/output  routines  are  presented  in 
the  software  and  are  used  to  read  and  store  images  in  conjunction  with  other 
image  processing  software,  digitizing  cameras,  and  output  display  devices. 

A design  is  never  finished  but  simply  stopped  due  to  time  limitations. 
Therefore,  this  software  does  not  comprise  a complete  set  of  image  processing 
functions  and  display  tools,  but,  rather,  provides  a good  base  for  most  appli- 
cations and  can  be  expanded  with  additional  future  image  processing  algorithms. 

The  rest  of  this  section  describes  the  basic  two-dimensional  digitized 
image  used  in  image  processing.  Section  2 describes  in  depth  the  image  pro- 
cessing functions  and  display  tools  contained  in  the  software.  Section  3.1 
describes  an  application  of  the  software  for  a hybrid  computer-opti cal  image 
recognition  experiment  (a  copy  of  the  full  paper  is  in  appendix  C).  Section 

3.2  describes  a preliminary  experiment  that  relates  the  near-field  and  far- 
field  scans  from  a single-mode  optical  fiber.  Section  4 contains  concluding 
remarks  and  future  applications.  Appendices  A and  B show  the  manual  and  source 
code  listing  for  the  software. 
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Figure  1-1.  American  Gothic  by  Grant  Wood  (1930),  oil  on  beaver  board 

(29- 7 / 8 11  x 251').  Friends  of  the  American  Art  Collection.  Copy- 
right, The  Art  Institute  of  Chicago.  All  rights  reserved. 

Figure  1-1  is  an  image  used  throughout  section  2.  This  painting  was 
chosen  because  it  is  a familiar  image  with  good  range  of  contrast,  shapes,  and 
detai 1 . 

This  report  is  based  entirely  on  the  Master  of  Science  thesis  of  Matthew 
B.  Weppner  (Electrical  and  Computer  Engineering  Department,  University  of 
Colorado,  Boulder,  1986).  It  has  been  edited  lightly,  and  the  format  has  been 
changed,  but  it  is  otherwise  intact. 

1.2  Image  Representation  and  Format 

A digitized  image  is  made  up  of  picture  elements  called  pixels.  A pixel 
is  a small  region,  usually  square,  whose  size  determines  the  resolution  of  the 
picture.  Each  pixel  is  assigned  a numerical  value  which  represents  the  aver- 
age light  intensity  over  the  area  of  the  pixel.  The  intensity  values  are  not 
absolute  measurements  but  represent  a relative  intensity  level.  In  this  soft- 
ware, each  pixel  value  must  be  a positive  integer,  one  byte,  which  allows  for 
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256  intensity  levels.  Possible  pixel  values  are  0 to  255,  and  the  maximum 
dynamic  range  over  the  256  intensity  levels  is 


maximum  dynamic  range  = 10  • log10  (— — ) = 24  dB.  (1) 

Although  this  represents  the  maximum  dynamic  range  of  the  software,  the 
dynamic  range  of  a digitized  image  is  determined  by  the  weakest  link  among  the 
input  digitizing  device,  software,  and  output  device. 

Each  digitized  two-dimensional  image  is  made  up  of  a finite  number  of 
pixels  with  a horizontal  width  and  vertical  height  (resolution  format)  speci- 
fied by  the  software.  The  image  is  represented  by  a two-dimensional  array, 
p(x,y),  which  is  shown  in  figure  1-2.  Each  element  p(x0,yQ)  represents  the 
intensity  value  of  a pixel  at  a specified  position  (x0,y0).  In  this  software 
valid  resolution  formats  are  256  x 240,  128  x 120,  64  x 60,  32  x 30,  and 
16  x 15. 

Figure  1-3  is  a 256  x 240  digitized  image  of  figure  1-1.  A CCD  (charge 
coupled  device)  camera  with  resolution  480  x 384  is  used  to  obtain  this  image. 
Four  similar  pictures  have  been  averaged  together  to  reduce  the  effects  of 
video  camera  noise  (section  2.3). 

2.  Image  Processing  Functions 
2.1  Single-Pixel  Functions 

Single-pixel  functions  are  those  in  which  the  (intensity)  value  of  a 
pixel  is  altered  according  to  an  algorithm  that  depends  on  the  original  value 
of  that  pixel.  The  functions  have  the  general  form, 

^xo*^o^  ^si ngl e^xo’^o^ ’ ' ^ ' 

where  p(xQ,y0)  is  the  current  pixel  being  processed  in  the  picture  p(x,y). 

These  functions  scan  the  entire  image  and  process  each  pixel  indivi- 
dually. In  general  these  functions  are  relatively  quick  to  execute  and  do  not 
require  complex  computation. 


3 


Width  = M 


0 


<-  x -> 


M-l 


Height  = N 


0 

t 

y 


N-l 


Figure  1-2.  Diagram  of  a M x N digitized  image  p(x,y). 


Figure  1-3.  256x240  digitized  image  of  figure  1-1. 
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2.1.1  Biasing 


Biasing  adds  an  input  scalar  value  to  each  pixel  of  a digitized  image. 
The  input  value  must  be  an  integer  within  the  valid  range  (0  to  255).  The 
general  form  of  the  equation  is 


where  $ = input  bias  value  and  has  the  range  -255  < 8 < 255. 

The  purpose  of  this  biasing  function  is  to  shift  the  spectrum  of  intens- 
ity levels  by  a specified  amount.  This  changes  the  average  intensity  level  of 
the  picture  and  can  alter  the  range  of  visual  contrast  of  an  image.  Image 
information  loss  results  if  computed  values  exceed  the  intensity  range  and  are 
truncated. 

Figure  2-1  shows  the  result  of  biasing  figure  1-3  with  8 = 50.  Fig- 
ure 1-3  initially  had  a good  range  of  exposure.  In  figure  2-1  the  picture 
appears  much  brighter,  but  the  brightest  parts  of  the  picture  appear  over- 
exposed because  the  pixel  values  in  those  areas  exceed  255.  Figure  2-2  shows 
the  result  of  biasing  figure  1-3  with  8 = -50.  In  this  case  the  picture  is 
much  darker  and  much  of  the  dark  area  appears  underexposed  because  the  pixel 
values  in  those  areas  have  values  that  are  less  than  0. 

2.1.2  Contrast  Enhancement 

Contrast  enhancement  increases  the  range  of  contrast  of  the  intensity 
values  of  a picture.  Pixel  values  within  an  input  range  less  than  256,  are 
linearly  expanded  to  the  full  range  of  0 to  255.  The  input  range  is  deter- 
mined by  a lowest  and  highest  input  bounds.  The  general  form  of  the  equation 


0 


255 


P(x0»y0)  + 8 < 0 

0 < p(VV  + 6 < 255 
P(x0>y0)  + 6 > 255 


(3) 


is 


f 


enhance 


0 a < 0 

(p(x0»y0) »a»b)  = a 0 < a < 255  (4) 
255  a > 255 
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Figure  2-1.  Biased  image  of  figure  1-3  (=+50). 


Figure  2-2.  Biased  image  of  figure  1-3  (=-50). 
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pcc 

where  a = (Ilf)  ( p(x0 »y0 )-a) , 

a = input  lower  bound, 
b = input  upper  bound, 

0 < a < b < 255. 

Pixel  intensity  values  not  within  the  range  of  the  lowest  and  highest  input 
values  will  exceed  the  valid  range  after  enhancement  and  are  truncated  to  0 
and  255,  respectively. 

Contrast  enhancement  makes  the  picture  appear  brighter  and  brings  out 
hidden  detail.  This  function  is  also  useful  for  expanding  the  dynamic  range 
of  intensities  of  a picture  for  other  image  processing  functions. 

Figure  2-3  shows  an  enhanced  image  of  Figure  1-3  using  the  values  a = 49 
and  b = 227,  which  are  the  lowest  and  highest  pixel  values  in  figure  1-3.  We 
call  this  process  "auto-enhance"  in  the  software  because  the  image  is  auto- 
matically adjusted  to  the  maximum  possible  contrast.  Figure  2-4  shows  another 
enhanced  image  of  Figure  1-3  using  a=100  and  b= 1 75.  This  image  is  deliber- 
ately over-enhanced  to  bring  out  the  detail  within  range  between  a and  b. 

Pixel  information  outside  this  range  is  lost. 

2.1.3  Negative  Formation 

Forming  a negative  of  a digitized  picture  is  analogous  to  creating  nega- 
tive image  in  photography.  The  range  of  pixel  intensities  is  reversed,  with 
the  brightest  pixels  becoming  the  darkest  and  vice  versa.  The  form  of  the 
equation  is 


^negati ve^P^xo*^o^  ~ P^xo5^o^* 

The  purpose  of  this  function  is  to  reverse  the  contrast  of  a picture  and 
create  an  inverse  image.  Figure  2-5  shows  the  negative  image  of  figure  1-3. 
This  function  is  helpful  in  bringing  out  detail  and  analyzing  certain  images. 
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Figure  2-3.  Enhanced  image  of  figure  1-3  (a=49,  b=227) 


Figure  2-4.  Enhanced  image  of  figure  1-3  (a=100,  b=175). 
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Figure  2-5.  Negative  image  of  figure  1-3. 

2.1.4  Thresholding 

Thresholding  is  used  to  set  a range  of  pixels  to  a desired  value.  Upper 
and  lower  pixel  input  bounds  are  used  to  establish  a range  of  intensity.  The 
general  form  of  the  function  is 

fthreshold(p(xo,yo) >a’b’T) 

,t  a < P(xQ,y0)  < b 

(6) 

( p(x  ,y  ) otherwise 

o o 

where  a = input  lower  bound, 
b = input  upper  bound, 

0 < a < b < 255, 
t = input  threshold  set  value, 

0 < x < 255. 

This  function  sets  a range  of  pixel  intensity  values  equal  to  an  input  thresh- 
old value.  This  can  be  used  to  distinguish  ranges  of  intensity  within  a pic- 
ture. 
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Figure  2-6  shows  the  thresholded  image  of  figure  1-3  with  the  values 
a = 0,  b = 128,  and  t = 64.  Contrast  information  within  this  range,  0 to  128, 
is  lost  and  diverts  our  attention  to  the  unthresholded  pixel  values. 

Figure  2-7  shows  an  image  of  figure  1-3  with  the  number  of  possible  pixel 
intensity  values  reduced  from  256  to  16.  This  is  done  by  setting  ranges  of 
pixels  into  16  threshold  levels,  by  thresholding  16  times  using  a = 0,  b = 15, 
and  t = 0,  then  a = 16,  b = 31,  and  t = 16,  etc. 

2.2  Multiple-Pixel  Functions 

Multiple-pixel  functions  are  image  processing  functions  which  operate  on 
one  pixel  as  a function  of  many  pixels  or  even  the  entire  digitized  image. 

The  functions  have  the  general  form. 


where  p'(x0,y0)  is  the  current  pixel  being  processed  in  the  picture  p(x,y). 

These  functions  scan  the  entire  picture  but  process  each  pixel  individ- 
ually. Since  each  pixel  is  a function  of  many  pixels  (or  the  entire  image)  we 
must  put  the  processed  results  into  a second  array,  p'(x,y).  This  is  neces- 
sary because  the  computer  is  capable  of  processing  only  one  pixel  at  a time, 
and  the  entire  original  image  is  needed  until  all  the  pixels  are  processed. 
Usually,  these  functions  are  relatively  slow  to  execute  and  require  more  com- 
plex computation  compared  to  single  pixel  functions. 

2.2.1  Convolution 

The  purpose  of  the  convolution  function  is  to  apply  a spatial  frequency 
filter  to  the  digitized  picture.  A straightforward  method  [1]  of  convolving 
two  functions  is  to  multiply  the  Fourier  transform  of  the  digitized  image  with 
a frequency  filter  function  in  the  frequency  plane.  The  inverse  Fourier 
transform  is  then  taken  to  transform  the  result  back  into  the  image  plane. 

The  general  form  is 


p ' (x  ,y  ) = f 
' ' o o'  i 


multiple^  o,Jo 


(x  ,y  ,p(x,y) ) , 


(7) 


P'(x,y)  = F_1{P(fx,fy)  • K(fx,fy)>, 


(8) 
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Figure  2-6.  Thresholded  image  of  figure  1-3  (a=0,  b= 128,  t=64). 


Figure  2-7.  Thresholded  image  of  figure  1-3  with  the  number  of  pixel 
intensity  values  reduced  from  256  to  16. 
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where  F{g(x,y)}  = G(f  ,f  ) = (forward)  Fourier  transform, 

x y 

F_1{G(x,y)}  = g(x,y)  = inverse  Fourier  transform, 

P(f  ,f  ) = F{p(x,y)}  = Fourier  transform  of  the  picture, 
x y 

K(fx,fy)  = frequency  filter  function. 

Computing  the  two-dimensional  Fourier  transform  of  the  picture  and  the 
inverse  transform  of  the  filtered  result  is  usually  very  time  consuming.  We 
can  also  filter  the  picture  with  respect  to  spatial  frequency  by  performing 
the  convolution  directly  in  the  image  plane.  Using  the  convolution  theorem  we 
obtain  [1], 


F{k(x,y)  * p(x,y)}  = P(fx>fy)  • K(fx,fy), 

k(x,y)  * p(x,y)  - F'i{P(fx,fy)  • K(fx,fy)}  - p‘(x,y)>  (9) 

-foo  -00 

where  f * g = / / g(c,n)g(x-c,y-n)  d^dn, 

— 00  —00 

k(x,y)  = F_1{K(f  ,f  )}  = convolution  kernel  function, 
x y 

and  * denotes  convolution.  Convolution  in  the  image  plane  is  the  mathematical 
equivalent  to  multiplication  in  the  Fourier  plane.  We  use  convolution  because 
it  is  simpler  to  implement  and  faster  when  the  kernels  are  small. 


In  image  processing  we  deal  with  digitized  pictures  in  matrix  form. 
Similarly,  we  must  chose  a finite  matrix,  k(x,y),  which  serves  as  the  convolu- 
tion kernel.  This  convolution  kernel  is  often  an  approximation  to  the  inverse 
Fourier  transform  of  the  frequency  filter  function,  K(fx,fy).  The  general 
form  of  the  two-dimensional  computer  convolution  formula  is 


where 


fconvolve(VVp(x’y)’k(x'y)) 


0 

a < 

0 

255 

a > 

255 

a 

otherwi se 

a = 0 + 


Nk 

( l 

n=l 


l k(m,n) 
m=l 


p(x0+m-mc,y0+n-nc)). 


(10) 


8 = bias  value, 

t)  = normalization  value, 
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l l k(m,n)  = 0 
n=l  m=l 

\ l l k(m,n)  otherwise 
n=l  m=l 

k(x,y)  = Mk  x Nk  matrix  (convolution  kernel, 

(mc,nc)  = center  of  matrix  k(x,y), 
mc  = (Mk  div  2)+l, 
nc  = (Nk  div  2)+l, 

and  "div"  denotes  integer  division,  which  ignores  the  remainder. 

Figure  2-8  shows  a low-pass  image  of  figure  1-3.  The  convolution  kernel 


k 


low-pass 


1 1 1 

111. 

1 1 1 


(11) 


This  kernel  is  a rough  approximation  to  a small  circles  in  the  Fourier  plane. 
The  result  is  loss  of  high  frequency  information  or  blurring. 

Figure  2-9  shows  a high-pass  filtered  image  of  figure  1-3.  The  convolu- 
tion kernel  is 


k 


high-pass 


0-10 

-14-1. 

0-10 


(12) 


This  kernel  is  a rough  approximation  to  a sinc(r)  function,  which  represents 
an  opaque  stop  in  the  Fourier  plane.  The  result  is  loss  of  low  frequency 
information  and  gives  an  image  with  only  the  edges  and  noise  remaining. 

2.2.2  Fourier  Transforms 

The  Fourier  transform  functions  transform  a digitized  image  in  either  the 
forward  or  reverse  direction.  This  is  similar  to  the  Fourier  transform! ng 
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Figure  2-8.  Low-pass  filtered  image  of  figure  1-3. 


Figure  2-9.  High-pass  filtered  image  of  figure  1-3. 
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properties  of  a lens,  although  a lens  performs  only  a forward  transformation. 
The  form  of  the  transform  functions  is 

P (^x >fy)  = fforward-FFT(P(x ) > 

P (x,y)  = f-j  nverse-FFT(p  (px  ,py)  ’^,e) » 

where  P(fx,fy)  = the  2-D  Fourier  transform  of  p(x,y), 
x,y  = spatial  variables  in  the  image  plane, 
f x , f y = frequency  variables  in  the  Fourier  plane 
8 = transform  plane  magnification  factor, 
e = transform  plane  scaling  factor. 

In  this  software  a digitized  image  is  assumed  to  be  a phaseless  intensity 
distribution.  To  process  the  Fourier  transform  we  convert  the  intensity 
distribution  into  a complex  amplitude  distribution.  For  the  forward  trans- 
form, 


(13) 

(14) 


P(VV  lacomplex^xo’^o^ 

^areal  ^xo *^o^  + ^aimaginary^xo,,yo^ 


(15) 


where  acomplex(x,y) 
Arbitrarily,  we  set 


= complex  amplitude  distribution  in  the  image  plane. 


areal(Vyo)  = <P(xo>yo))  k' 


a.  . (x  ,y  ) = 0. 
imaginary'  o o' 


(16) 


Similarly  for  the  inverse  transform, 


P(fx  .fx  ) = |A 

xo  xo 


x_ ’ x_'  i complex'  x 5 y 


(f.  .fy 

0 Jo 


= (A  , (f  ,f  ))2  + (A.  . (f  ,f  ))2, 

' real'  xQ*  yQ " ' imaginary'  xQ*  y0" 


(17) 


where  Acomp]ex(fx,fy)  = complex  amplitude  distribution  in  the  Fourier  plane. 


Again,  we  set 
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(18) 


Areal<fx  -fy  > = <P<fx  >fy  » 
0 J0  0 Jo 


72 


^imaginary^x^y^ 


To  describe  the  transforming  process  in  detail  we  start  with  the  basic 
Fourier  transform  equations.  The  equation  for  the  forward  Fourier  transform 
is  [1] 

+°0  +00  — j 2tT  (f  X+f  V ) 

F{g(x,y)}  = / / g(x,y)  • e r dxdy  = G(f  ,f  ),  (19) 

— 00  — OO  J 

and  the  inverse  Fourier  transform  is 

+oo  +oo  +j2ir(f  x+f  y) 

F'*{G(f  f )}  - / / G(f  f ) • e x y dfdf  -g(x,y).  (20) 

* — 00  — 00  * * 


Since  we  are  interested  in  performing  the  two-dimensional  Fourier  transform  on 
a digitized  picture,  we  use  the  Fast  Fourier  Transform  algorithm  (FFT)  to  take 
advantage  of  its  inherent  speed,  which  is  necessary  for  the  transformation  of 
the  relatively  large  digitized  pictures.  The  FFT  algorithm  in  its  basic  form 
transforms  only  a one-dimensional  array.  To  perform  the  Fourier  transform  in 
two  dimensions  we  rewrite  the  two-dimensional  transform  as  successive  one- 
dimensional transforms: 


+°°  +°°  -2irf  x — j 27r f y 

Ffg(x,y)}  = f f g(x,y)  • e x • e y dxdy 

+oo  —00 

+oo  +oo  -2irf  x — 2tt f y 

= / ( / g(x,y)  • e x dx)  • e y dy 

— 00  -00 

= Fy{Fx{g(x,y) }}  = Fx(Fy{g(x,y)}}. 
Similarly  for  the  inverse  transform. 


F"1fg(fx.Fy»  = Fy-i{Fx'i{g(fx)fym 


Fx'ltF/>tg(fx.fy))}. 


(21) 


(22) 
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Thus,  to  obtain  the  two-dimensional  transforms  we  can  first  transform  all  the 
(vertical)  columns,  and  then  transform  all  the  (horizontal)  rows,  or  vice 
versa. 

If  the  input  picture  function  has  moderate  size  relative  to  the  bounds  of 
the  picture  array,  the  transformed  result  will  occupy  a small  portion  of  the 
transform  plane  and  will  not  contain  many  nonzero  data  points  or  pixels.  To 
achieve  a transformed  result  of  moderate  size  with  an  input  of  moderate  size 
we  use  transform  plane  magnification.  This  results  in  a transform  with  more 
nonzero  elements,  though  there  is,  in  reality,  no  more  information  than 
before.  Its  usefulness  is  mainly  for  visual  display. 

Transform  plane  magnification  involves  shrinking  the  input  plane  to 
achieve  a stretching  or  magnification  in  the  transformed  plane.  Using  the 
similarity  theorem  for  Fourier  transforms  we  see  that  [1] 

1 fx  fv 

F{g(ax,by)}  = -R-rrr  • G b )•  (23) 

Thus,  for  the  magnified  forward  and  inverse  transforms, 

Ff  p(f,f)}  = s2  • P(sfx,efy),  (24) 

f f 

F'M  - e2  • P(sx,sy),  (25) 


where  8 = transform  plane  magnification  factor.  To  obtain  a shrinking  of  the 
input  plane  without  loss  of  information,  the  input  plane  is  enlarged  by  a 
factor  of  8 and  padded  with  0's.  This  creates  a factor-of-8  magnification  in 
the  transform  plane  with  no  loss  of  information  in  the  input  plane  because  the 
input  is  only  effectively  shrunk  with  the  use  of  larger  FFT  arrays.  The 
penalty  is  increased  computing  time. 

One  of  the  requirements  of  the  FFT  algorithm  is  that  the  number  of  ele- 
ments in  the  one-dimensional  input  array  be  a power  of  2.  Therefore,  for  nor- 
mal situations  (8  = 1),  the  number  of  elements  in  each  one-dimensional  FIT 
array  must  also  be  a power  of  2.  In  both  the  horizontal  and  vertical 


17 


directions,  the  size  of  the  FFT  array  must  equal  the  number  of  pixels  in  that 
direction.  If  the  number  of  pixels  is  not  a power  of  2,  the  number  of  ele- 
ments used  for  the  FFT  array  must  be  the  next  higher  power  of  2.  For  an  M x N 
pi cture, 


a 


max 


flog 

intMog  V 


for  horizontal  transforms. 


max 


rl23Jh 

intMog  2' 


for  vertical  transforms, 


(26) 


where  2a  = minimum  number  of  elements  in  the  FFT  array,  and 
max-jnt()  = rounds  up  to  the  next  integer. 

To  adjust  for  the  Fourier  plane  magnification,  we  use 


A = 8 • 2a 


(27) 


where  A = number  of  elements  in  the  FFT  array  to  be  used,  and 
8 = 2n  for  n > 0. 

Since  a is  an  integer,  8 must  be  a power  of  2 to  ensure  that  A is  also  a 
power  of  2. 

Since  the  number  of  elements  used  in  either  direction  usually  exceeds  the 
numbers  of  pixels,  extra  FFT  array  elements  are  set  to  0 before  being  trans- 
formed. After  transformation,  the  extra  FFT  array  elements  are  ignored. 

We  are  careful  in  choosing  an  input  image  that  contains  0's  around  the 
boundaries  of  the  two-dimensional  picture  array  to  prevent  ringing  due  to  dis- 
continuous edges.  Also,  a convention  set  by  the  FFT  algorithm  requires  that 
the  zero  frequency  term  (the  dc  term)  be  shifted  to  the  first  element  of  the 
array  to  be  transformed.  The  transformed  result  also  assumes  that  the  first 
element  is  the  dc  term.  Therefore,  the  horizontal  rows  and  vertical  columns 
are  shifted  to  move  the  center  dc  term  from  the  assumed  center  of  the  picture 
(appendix  A)  to  the  upper  left  position  (0,0).  After  the  transformati on  the 
rows  and  columns  are  shifted  once  again  to  move  the  dc  term  back  to  the  center 
position. 
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Once  the  complex  amplitude  distribution  is  obtained  the  array  is  trans- 
formed. The  vertical  columns  of  the  picture  are  transformed  first  and  the 
horizontal  rows  are  transformed  second.  For  the  forward  transform, 

Acomplex(fx’fy)  = FFTx{FFTy{acompiex(x,y)}} , (28) 

and  for  the  inverse  transform, 


acomplex(x»y)  FFTx  ^FFTy  1{Acompiex(fx,fy)}}  • (29) 

Now  the  transformed  amplitude  distribution  must  be  converted  back  into  an 
intensity  distribution.  For  the  forward  transform, 


P^xQ,fy0)  (Veal  ( V0» *y0^2+(  Vmagi nary(V0» V0^2  * 


(30) 


and  for  the  inverse  transform. 


P(xo»y0)  (areal  (xo’^o^2  + (aimaginary(xo’^o^2 * (^1) 

Since  phase  is  not  retained  after  converting  back  to  an  intensity  distri- 
bution, these  transformed  images  are  not  completely  analogous  to  amplitude 
transforms,  which  have  phase.  Consequently,  the  Fourier  transform  functions 
described  here  are  not  conservative,  and  we  cannot  undo  the  transform  process. 
Using  the  Fourier  integral  theorem  for  amplitude  distributions  [1]  we  see  that 

g(x,y)  = F" 1 {F{g (x ,y ) } } = F {F “ 1 (g(x,y)} . (32) 

But,  since  we  are  assuming  intensity  distributions  with  no  retention  of  phase 
terms. 


P(x,y)  * fj nverse-FFT^forward-FFl(P(x )) ) 

* Vorward-FFT^  i nverse-FFl(P(x,y) ) ) * (33) 

We  must  be  careful  in  interpreting  these  transform  functions,  keeping  in  mind 
that  the  digitized  input  pictures  and  transformed  results  are  intensity  dis- 
tributions and  have  certain  limitations. 
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Since  Fourier  transforms  usually  have  large  peaks  we  use  an  integer  scal- 
ing factor,  e,  to  enhance  smaller  details  which  would  not  normally  seen 
because  of  the  limited  dynamic  range  of  the  pixel  values  (0  to  255).  This 
scaling  process  is  similar  to  overexposing  film  to  observe  detail  normally 
drowned  out  by  the  strong  peaks.  The  scaling  factor,  e,  represents  the  level 
of  overexposure  to  which  the  transformed  values  are  assigned  pixel  values.  A 
scaling  factor  of  1 means  the  highest  transformed  FFT  value  is  assigned  255 
and  the  rest  of  the  values  are  scaled  accordingly.  For  values  of  e greater 
than  1 the  scale  is  simply  multiplied  by  that  factor.  Pixel  values  which  are 
computed  to  be  higher  than  255  are  truncated.  For  the  forward  transform, 

(p  0 < p < 255 
P(f  ,f  ) = ] 

xo  yo  (255  otherwise 

P = • P(fx  .fy  ).  (34) 

" xo  yo 

where  e = integer  scaling  factor  (e  > 1), 

h = highest  value  contained  in  P(fx,fy). 


For  the  inverse  transform. 


p(x  ,y  ) = 
o o 


p = ( 


(255 
255  • s 


0 < p < 255 
otherwise 

P(x0.y0). 


(35) 


where  e = integer  scaling  factor  (e  > 1), 

h = highest  value  contained  in  p(x,y). 


Figure  2-10  shows  a 128  x 120  Fourier  transform  of  figure  1-3.  128  x 120 
is  the  highest  resolution  that  can  be  transformed  by  our  computer,  since  this 
resolution  requires  120  kbytes  of  memory  for  the  image  alone  using  complex 
floating  point  arithmetic  which  requires  8 bytes  per  pixel  (a  256  x 240  image 
requires  480  kbytes  of  memory).  Since  the  image  of  figure  1-3  is  large 
relative  to  the  bounds  of  the  image  array,  the  resulting  transform  is  quite 
small.  Enlargement  is  achieved  with  a Fourier  plane  magnification  factor 
8=4.  To  bring  out  the  finer  transform  detail  the  image  was  scaled  or  over- 
exposed by  a factor  of  e = 3000. 
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Figure  2-10.  Fourier  transform  of  figure  1-3. 

Figure  2-11  shows  a 128  x 120  image  of  a circ(r)  function  with  r = 20  and 
height  = 255.  Figure  2-12  shows  the  Fourier  transform  of  figure  2-11,  which 
has  the  form  sinc2(r).  Once  again  a Fourier  plane  magnification  factor  B = 4 
is  used.  To  bring  out  the  detail  in  the  secondary  maxima  a scaling  factor 
e = 50  is  used. 

2.2.3  Magnification/Demagnification 

Magnification  increases  the  physical  scaling  of  a digitized  picture  by  a 
factor  of  2 in  both  the  horizontal  and  vertical  directions.  A factor  of  2 is 
the  only  value  provided  by  the  software  to  ensure  exact  processing  without 
aliasing  between  pixels.  This  process  is  similar  to  2X  magnification  using 
lenses.  Information  contained  in  the  picture  is  decreased  by  a factor  of  4 
because  only  the  center  fourth  of  the  area  of  the  image  is  capable  of  magnifi- 
cation. Information  outside  the  area  to  be  magnified  is  lost.  Magnified 
pixels  are  simply  replicated  to  produce  the  larger  image.  The  general  form  of 
the  equation  is 
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Figure  2-11.  128x120  computer  generated  circ(r)  function  with  r=20. 


Figure  2-12.  Fourier  transform  of  figure  2-11. 
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(36) 


fmagnify(xo>VP(x-'>')>  = P<xa+(xo  div  d1v  2>> 

where  xa  = M div  4, 
ya  = N div  4, 

M x N = picture  resolution  format. 


Figure  2-13  shows  the  magnified  image  of  figure  1-3. 


Demagnification  decreases  the  scaling  of  a digitized  picture  by  a factor 
of  2 in  both  directions.  Information  in  the  picture  is  decreased  by  a factor 
of  4 because  averaging  of  four  pixels  is  used  to  create  each  demagnified 
pixel.  This  is  similar  to  2X  image  reduction  using  lenses.  Information  taken 
from  beyond  the  bounds  of  the  input  image  is  obviously  unavailable  and  the 
demagnified  pixels  are  set  to  0.  The  general  form  of  the  equation  is 


^demagnify  (xo’^o’^x’^  j g 


< x < x.  , y < y < y, 
o b Ja  Jo 


otherwise 


1 2 2 

« - i i x 

j=i  i-i 


P(2(x0  - xa)  + 1 ,2(y0  - ya)  + j). 


(37) 


where  xa  = M div  4 x^  = 3(M  div  4 ) “ 1 , 

ya  = N div  4 y^  = 3(N  div  4)"1, 

M x N = picture  resolution  format. 


Figure  2-14  shows  the  demagnified  image  of  figure  1-3.  Both  functions  can  be 
used  repeatedly  on  an  image  to  achieve  higher  magnification  or  demagnification 
factors. 


2.2.4  Noise  Filtering 

Noise  filtering  is  used  to  reduce  random  noise  such  as  video  camera  noise 
or  noise  generated  by  other  image  processing  functions.  The  algorithm  for  the 
function  is  from  reference  [2].  When  the  intensity  value  of  a pixel  differs 
from  the  average  of  its  eight  surrounding  neighborhood  pixels  by  more  than  a 
set  value,  e,  the  pixels  are  set  to  the  average  neighborhood  value.  The 
general  form  of  the  function  is 
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Figure  2-13.  2X  magnified  image  of  figure  1-3. 


Figure  2-14.  2X  demagnified  image  of  figure  1-3. 
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otherwise 


(38) 


where  e = input  set  value, 

Plt.  . .,P8  = the  eight  neighboring  pixels  of  p(xQ,y0). 

The  average  of  the  neighboring  pixels  of  this  function  is  obtained  by  automa- 
tically convolving,  in  software,  the  image  with  the  noise  filtering 
convol ution  kernel , 


The  difference  between  p(xQ,y0)  and  the  average  is  then  compared  with  e,  and 
p(x0,y0)  is  set  accordingly  by  the  function. 

Figure  2-15  is  an  image  of  figure  1-3  with  randomly  generated  pixels  of 
value  255  added.  This  random  noise  is  similar  to  some  types  of  video  camera 
noise  under  extreme  conditions.  Figure  2-16  is  an  image  of  figure  2-15  with 
the  noise  filtering.  Looking  closely  at  the  cleaned  image  we  can  see  that  the 
noise  is  greatly  reduced  but  not  completely  eliminated. 

2.2.5  Resolution  Changing 

The  resolution  changing  function  changes  the  hori zontal -by-vertical  reso- 
lution format  of  the  input  image.  Changing  resolution  does  not  alter  the 
scale  of  the  picture  but  alters  only  the  number  of  pixels  in  both  the  horizon- 
tal and  vertical  directions. 

Decreasing  resolution  decreases  the  number  of  pixels  by  a factor  of  4n  by 
reducing  the  number  of  pixels  in  both  directions  by  a factor  of  2n.  This 
results  in  a loss  of  information  because  groups  of  4n  pixels  are  averaged  to 
produce  pixels  in  the  lower  resolution  picture.  The  general  form  of  the  equa- 
tion is 


(39) 


Figure  2-15.  Image  of  random  noise  superimposed  on  figure  1-3. 


Figure  2-16.  Reduced  noise  image  of  figure  2-15. 
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(40) 


^dec  rea  se  - res  ^o’^o’P^5^ ) ,n) 

2n 

= \ l l P(2n  • x0  + 1 >2n  * yQ  + J)  > 

4n  j=l  i=l  0 0 

where  n = resolution  reduction  factor  (n  s>  1).  The  new  resolution  format  is 

M' =^andN'=7*  <«> 

where  M x N = resolution  format  of  p(x,y), 

M'x  N'=  resolution  format  of  p'(x,y). 

Figure  1-3  is  the  highest  resolution  format  handled  by  the  software  at 
256  x 240.  Figures  2-17,  2-18,  2-19,  and  2-20  are  reduced  resolution  images 
of  figure  1-3  with  respective  resolution  formats,  128  x 120,  64  x 60,  32  x 30, 
and  16  x 15. 

Increasing  resolution  increases  the  number  of  pixels  by  a factor  of  4n  by 
increasing  the  number  of  pixels  in  both  directions  by  a factor  of  2n.  This 
results  in  no  gain  of  information  because  pixels  are  simply  replicated  by  a 
factor  of  4n.  The  general  form  of  the  equation  is 

fincrease-res(xo»yo’P(x’y)’n)  = P(xo  div  ^o  div  2")’  (42) 

where  n = resolution  enlargement  factor  (n  > 1).  The  new  resolution  format  is 

M'  = 2n  • M and  N'  = 2n  . N,  (43) 

where  M x N = resolution  format  of  p(x,y), 

M'x  N'=  resolution  format  of  p'(x,y). 

Increasing  resolution  does  not  change  appearance  of  an  image  but  may  be 
necessary  for  matching  resolution  formats  of  different  imaging  software  or 
input/output  devices. 
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Figure  2-17.  128x120  image  of  figure  1-3. 


Figure  2-18.  64x60  image  of  figure  1-3. 
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Figure  2-19. 


32x30  image  of  figure  1-3. 


Figure  2-20. 


16x15  image  of  figure  1-3. 
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2.3  Two-Picture  Functions 


Two-picture  functions  are  image  processing  functions  that  compute  each 
pixel  as  a function  of  the  corresponding  pixels  of  two  pictures,  p^x.y)  and 
p (x,y) . The  functions  have  the  form, 


Pi(xo^o)  ^two-picture^i  (xo’^o^  *^2  ^xo5^o^  » 

where  P^Xq^q)  is  the  current  pixel  being  processed  in  the  picture  p (x,y). 
(In  this  algorithm,  p (x,y)  is  replaced  by  the  new  array.)  Examples  of  two- 
picture  functions  include  averaging,  addition,  subtraction,  multiplication, 
and  logical  functions.  The  averaging  function  is  the  only  two-picture  func- 
tion included  in  the  software  and  is  a very  useful  method  of  noise  reduction. 

These  functions  scan  the  entire  picture  and  process  each  pixel  individ- 
ually. Both  pictures  must  have  the  same  resolution  format.  In  general  these 
functions  are  quick  and  usually  do  not  require  complex  computation. 


2.3.1  Averaging 


Averaging  is  creating  a picture  that  is  the  pixel -by-pixel  average  of  two 
pictures.  The  general  form  of  the  equation  is 


^ average^1  ^xo’^o^  *P2 (xo*^o^ 


MVV  + P2(x0,y0) 


(45) 


where  p (x,y)  = current  picture  in  memory, 
p (x,y)  = input  picture  file. 


The  purpose  of  this  process  is  to  reduce  random  noise  such  as  video  camera 
noise.  Averaging  two  or  more  pictures  can  reduce  this  type  of  noise 
si gnificantly. 

Using  only  a two-picture  function,  we  may  also  average  any  number  of  pic- 
tures that  is  a power  of  2.  For  example,  four  pictures  can  be  averaged  in 
three  steps  as  follows. 
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P3(xo’y0)  ~ ^average^a ^xo’^o^ ’^4 ^xo’^o^  5 
Pi^xo’^o^  ” ^average^i  ^xo’^o^  * ^2  ^xo’^o^  ’ 

Pi(xo^o)  ” ^average^i  ^xo’^o^  *P3  (xo’^o^  * 

where  p^x.y)  becomes  the  average  of  the  four  pictures  p (x,y),  p (x,y), 
P3(x,y),  and  p4(x,y). 

Figure  1-3  shows  the  result  of  averaging  four  images  similar  to  fig- 
ure 2-21.  The  eye  is  a marvelous  averaging  device;  the  image  looks  fine 
viewed  live  on  a video  screen  but  is  noisy  when  digitized  and  frozen,  as  in 
figure  2-21.  Figure  2-22  shows  a noisy  image  of  the  output  of  a single-mode 
optical  fiber.  The  image  is  obscured  by  randomly  generated  noise  as  a result 
of  the  low-light  conditions.  Figure  2-23  is  the  average  of  eight  images  taken 
separately  but  all  similar  to  figure  2-22.  The  averaged  image  contains  much 
less  random  noise. 


Figure  2-21.  Unaveraged  image  of  figure  1-3. 
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Figure  2-22. 


Noisy  image  of  a Gaussian  beam. 


Figure  2-23. 


Averaged  image  of  figure  2-22. 
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2.4  Image  Information  Display  Tools 

Image  Information  Display  Tools  compute  and  display  information  about  a 
digitized  image.  Some  of  these  tools  are  custom  in  that  they  are  designed  to 
meet  the  needs  of  a specific  user.  The  Gaussian  and  g-profile  curve-fitting 
tools  are  specially  designed  to  meet  the  needs  of  characteri zi ng  beam  spot 
radius  of  single-mode  and  multimode  optical  fibers.  The  hi stogram-of- 
intensity  graph,  line  intensity  graph,  and  three-dimensional  intensity  plot 
have  much  more  universal  applications. 

2.4.1  Gaussian  Curve  Fitting 

Two-dimensional  Gaussian  intensity  distributions  are  found  in  lasers, 
optical  fibers,  and  integrated  optical  devices.  The  purpose  of  the  Gaussian 
curve  fitting  tool  is  to  take  a one-dimensional  line  of  data  from  a digitized 
image  and  fit  a Gaussian  curve  to  it. 

Since  the  dark  pixel  value  of  most  video  cameras  is  not  0,  we  usually 
want  to  bias  the  image  with  the  negative  of  the  dark  level,  that  is,  to  clip 
the  dc  level.  This  dc  level  may  be  taken  to  be  the  peak  value  of  the  intens- 
ity histogram.  Also,  since  these  images  are  usually  taken  under  low-light 
conditions  we  fit  to  the  Gaussian  function  using  only  pixel  values  greater 
than  a noise  cutoff  value,  n,  which  is  taken  to  be  the  maximum  absolute  value 
of  the  video  and  quantization  noise.  We  subjectively  set  the  value  of  n by 
looking  at  a line  intensity  scan  of  the  beam  image.  After  removal  of  the  dark 
level,  pixel  values  less  than  or  equal  to  n are  ignored  in  the  curve  fitting 
routine. 

After  biasing  and  deciding  on  an  appropriate  n,  we  calculate  the  centroid 
or  first  moment  of  the  image.  The  centroid  is  the  center  position  of  a digit- 
ized picture  whose  coordinates  are  determined  by  the  weighted  average  along 
each  axis.  The  general  form  of  the  coordinates  of  the  centroid  [3]  is 

N M N M -1 

c = ( l I x • p ' (x,y)  ) • ( l 7 p 1 (x,y)  ) , 

y=0  x=0  y=0  x=0 

N M N M -1 

C = ( I l y • p'(x,y)  ) • ( £ l p'(x,y)  ) , (47) 

3 y=0  x=0  y=0  x=0 
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where  (cx,Cy)  = centroid  of  image, 

(p(x,y)  p(x,y)  > n 
P'(x,y)  = l 

to  otherwise, 

n = noise  cutoff  value. 

The  general  form  of  the  one-dimensional  Gaussian  curve  is 

P(u)  - A • e-a<u-c)\ 


where  u,c  = 


x, c  for  horizontal  line  fit  at  y = c 

x y 

y, c  for  vertical  line  fit  at  x = c . 

y x 


(48) 


To  fit  a Gaussian  curve  to  the  data,  we  linearize  by  taking  logarithms  and  use 
a linear  least  squares  fit.  In  linear  form, 

ln(p(u) ) = ln(A)  - ct(u-c)2. 


y 1 = 1 n(A)  - ax' , 


(49) 


where  y'  = ln(p(u))  and  x'  = (u-c)2.  We  now  have  the  basic  equation  of  a 
line. 


y'  = mx'  + b (50) 

where  m = -a  and  b = ln(A).  We  use  the  method  of  least  squares  to  determine 
the  line  of  best  fit  with  the  formulas  [4], 

k l x1  ,yi ’ - Jxj'  l y,.' 
m = ~ k l - (I  x,  ‘)2  > 

l (x,')2  l yf'  - l xi'  l xi'yi, 
b = k l (x^)2  - (J  x^)2  “ (51) 


where  y^ 1 = 1 n(p(u) ) 
x^ ' = (u-c)2 


but  not  used  if  p(u)  < n, 
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k = number  of  data  points  for  which  p(u)  > n 


< M for  horizontal  fit 


< N for  vertical  fit. 


After  the  fit  is  completed  we  have 


A = eb  and  a = -m. 


(52) 


We  use  these  parameters  to  calculate  the  beam  spot  radius  of  the  fitted 
data.  The  beam  spot  radius,  r,  is  defined  as  the  1/e2  point  relative  to  the 
maximum  of  the  intensity  distribution.  Using  the  a parameter,  we  see  that 


Examples  of  the  Gaussian  curve  fitting  tool  can  be  seen  in  section  3.2. 

2.4.2  Multimode  Curve  Fitting  (g-Profile) 

The  purpose  of  the  multimode  curve  fitting  tool  is  to  fit  a power-law 
profile,  known  as  a g-profile,  to  a graded  index  multimode  fiber  with  uniform 
cladding.  The  index  profile  is  given  as  [5] 


where  A = (nQ2  - n22)/2nQ2. 

This  equation  is  made  up  of  four  parameters,  nQ,  A,  core  radius  a,  and  expo- 
nent g.  The  parameters  a and  g are  of  most  importance. 

Since  fitting  a set  of  data  to  four  parameters  is  not  trivial,  a simpli- 
fication is  needed.  Following  Cherin  [6],  we  use  the  approximation, 


1/e2  = e-“r2 


and 


(53) 


n(r)  = n0[l  - 2i(r/a)g]1/2 


(54) 
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A n(r)  = An [1  - (r/a)g]. 


(55) 


where  An  = nL  - n2  and  An(r)  = n(r)  - n2. 

This  new  equation  now  contains  only  three  parameters.  An,  a,  and  g. 

Since  we  are  interested  in  looking  at  digitized  images  of  the  near  field 
of  multimode  fibers,  the  data  are  in  the  form  of  an  intensity  distribution, 

I(r)  = I0[l  - (r/a)9].  (56) 

Taking  logarithms  of  both  sides,  we  see  that 

1 n [1  - I ( r ) / I 0 ] = g • ln[r]  - g • ln[a],  (57) 

which  follows  the  basic  equation  of  a line, 

y = mx  + b,  (58) 


where  y = ln[l  - I(r)/I  0 ] , 
m = g, 
x = ln[r], 
b = -g  • 1 n [a] . 

Using  a least-squares  fit,  eq  (51),  we  can  solve  for  m and  b for  a given 
IQ.  The  g-parameter  we  seek  is  the  slope  m,  and  the  core  radius  is 


a = 


-b/m 

e 


(59) 


To  find  the  center  of  the  two-dimensional  profile,  a calculation  of  the 
centroid  is  performed,  eq  (47),  similar  to  the  Gaussian  curve  fit  in  section 
2.4.1.  From  the  coordinates  of  this  centroid  we  can  determine  the  center 
line,  in  either  horizontal  or  vertical  direction,  and  the  center  of  that  line 
from  which  we  compute  the  radius. 

In  the  EIA  (Electronic  Industries  Association)  standard  FOTP-58  [7],  only 
the  data  points  in  the  10  percent  and  80  percent  range  of  the  total  profile 
are  used.  In  other  words,  valid  data  points  used  in  the  curve  fitting  are 

only  those  which  fall  between  0 . 1 1 o anc*  0. 81 0 - 
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I cannot  be  obtained  by  a least-squares  fit  with  only  two  parameters. 
Consequently,  we  optimize  the  fit  with  respect  to  I by  looking  at  the  quality 
of  the  fit,  which  is  taken  as  the  error, 

E = l (I (r)  - I0[l  - (r/a)9])2.  (60) 

We  initially  set  I to  be  equal  to  Imax,  the  maximum  pixel  intensity  in  the 
digitized  image.  We  then  minimize  the  error  function,  E,  by  monitoring  E as 
I is  changed  and  recalculating  the  parameters  a and  g.  Once  E is  minimized, 
a final  set  of  parameters  for  I , a,  and  g is  reached. 

2.4.3  Intensity  Histogram 

An  intensity  histogram  is  a two-dimensional  graph  of  relative  frequency 
versus  the  pixel  intensity  values,  0 to  255.  The  relative  frequency  is 
defined  as  the  number  of  occurrences  of  a pixel  intensity  on  a normalized 
scale.  This  scale  is  defined  to  be  1.0  at  the  most  numerous  pixel (s)  and 
scaled  accordingly  for  the  remaining  pixel  values. 


I N I EH  S I T ¥ 

Figure  2-24.  Intensity  histogram  of  figure  1-3. 
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F'gure  2-25.  Line  '"tensity  scan  of  figure  1-3. 
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Figure  2-25.  _ine  intensity  scan  of  rigure  2-11. 
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Figure  2-27. 


Three-dimensiona1  intensity  plot  of  ~'g-re  1-3. 


Figure  2-28.  Three-dinensi  ona1  inte^sit)  pi  ct  :f  *'■; 
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Using  a histogram  we  can  look  at  the  relative  contrast  range  of  an  image. 
Also  computed  are  the  highest,  lowest,  average,  and  highest  relative  frequency 
(or  peak)  pixel.  Figure  2-24  shows  a histogram  of  figure  1-3.  The  highest 
pixel  value  is  227,  the  lowest  is  49,  the  average  is  121,  and  the  peak  is  79. 

2.4.4  Line  Intensity  Scan 

The  line  intensity  scan  tool  allows  us  to  graph  the  pixel  intensity 
values  across  a specified  line  in  either  the  horizontal  or  vertical  direction. 

Figure  2-25  shows  a horizontal  intensity  scan  of  line  120  of  figure  1-3. 
Figure  2-26  shows  a horizontal  intensity  scan  across  line  120  of  figure  2-11. 
Using  these  scans  we  can  analytically  monitor  the  results  of  image  processing 
functions  and  read  image  data  with  reasonable  accuracy. 

2.4.5  Three-Dimensional  Intensity  Plot 

The  three-dimensional  intensity  plot  generates  a visual  surface  whose 
height  is  determined  by  the  intensity  levels  of  an  image  versus  two- 
dimensional  positions  [8], 

Figure  2-27  shows  a three-dimensional  plot  of  figure  1-3.  Figure  2-28 
shows  a three-dimensional  plot  of  figure  2-11.  Using  these  plots  we  can 
better  visualize  the  image  intensity  distributions. 

3.  Applications  of  the  Software 

3.1  Creating  a Model  Image  for  a Hybrid  Computer-Optical 
Processing  Experiment 

The  purpose  of  our  experiment  in  hybrid  computer-optical  processing  was 
to  do  real-time  pattern  recognition  at  video  rates.  The  computer  was  respon- 
sible for  preprocessing  a model  image.  This  model  was  then  permanently  stored 
as  a Fourier  transform  hologram.  Optical  processing  was  then  used  in  real 
time  to  recognize  a live  image.  For  complex  enough  images,  optical  processing 
will  be  much  faster  than  the  computer. 
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In  this  experiment  we  created  a model  image  for  real-time  pattern  recog- 
nition [9].  A liquid-crystal  television  (LCTV)  was  used  as  a real-time  trans- 
parency to  display  digitized  images.  We  used  an  optical  processor  in  series 
with  a converging  beam  processor  [1]  to  display  the  Fourier  transform  of  a 
specially  prepared  model  and  generated  its  Fourier  transform  hologram.  Actual 
images  were  then  displayed  on  the  LCTV,  and  their  (optical)  Fourier  transforms 
illuminate  the  hologram.  The  convolution  theorem,  eq  (9),  shows  that  we  can 
obtain  the  convolution  of  the  model  with  the  actual  image  as  our  method  of 
image  recognition.  The  purpose  of  this  experiment  was  to  use  the  computer's 
advantage  in  digital  image  input/output  processing  and  also  take  advantage  of 
optical  processing's  superior  convolution  speed. 

As  a test  256x240  model  image,  we  used  an  ordinary  pair  of  pliers  as  in 
figure  3-1.  We  chose  the  pliers  for  their  distinct  shape  and  applicability  to 
image  recognition  in  manufacturing.  The  shiny  surface  of  the  pliers  required 
that  we  use  very  diffuse  and  uniform  illumination  with  a black  background  to 
bring  out  the  plier's  overall  shape  rather  than  accent  the  reflected  light 
from  its  contours. 

Since  the  background  of  figure  3-1  is  not  perfectly  black,  the  noise 
cleaning  function  of  section  2.2  was  used  to  remove  noise  and  spurious  bright 
spots.  To  remove  the  background  entirely,  we  used  the  thresholding  function 
and  specified  a threshold  range  to  set  that  area  to  0,  being  careful  not  to 
erase  any  detail  of  the  pliers.  The  resulting  image  is  shown  in  Figure  3-2. 

To  make  the  outline  shape  of  the  pliers  as  distinct  as  possible  we  used 
the  thresholding  function  once  again  and  set  all  the  pixels  of  the  pliers  to 
255.  The  resulting  image  is  shown  in  figure  3-3. 

From  the  model  image  we  wished  to  obtain  the  edge  pattern  of  figure  3-3. 
To  do  this  we  normally  convolve  the  object  with  a high-pass  convolution  ker- 
nel, given  in  eq  (12).  In  this  case  the  image  will  be  projected  on  an  LCTV 
with  resolution  of  only  140  x 120  pixels.  Since  this  lower  resolution  format 
does  not  match  exactly  any  of  the  formats  handled  by  this  software,  the  image 
resolution  of  the  model  image  remained  at  256  x 240  to  achieve  maximum  detail. 
To  compensate  for  the  larger  pixels  of  the  LCTV  we  used  a larger  convolution 
kernel , 
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Figure  3-1.  256x240  digitized  image  of  pliers. 


Figure  3-2.  Background  processed  image  of  figure  3-1. 
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model 


0 0 0 -1  0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

-1  0 0 4 0 0 -1 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 -1  0 0 0 . 


(61) 


This  convolution  kernel  is  analogous  to  eq  (12),  but  produces  an  outline  of 
the  image  that  is  three  times  wider,  which  is  large  enough  for  the  LCTV.  The 
final  model  image  is  shown  in  figure  3-4.  Details  of  the  pattern  recognition 
experiment  are  given  in  reference  [9]  and  are  reproduced  as  appendix  C. 


3.2  Analyzing  Near-  and  Far-Field  Images  for  an  Integrated 
Optics  Experiment 

The  propagation  of  light  in  the  form  of  a Gaussian  beam  has  great  inter- 
est in  the  fields  of  optical  fibers  and  integrated  optics.  The  important 
parameters  of  a Gaussian  beam  are  the  radius  of  curvature,  R(z),  of  the 
spherical  wave  and  its  spot  size,  w.  With  these  parameters  we  can  predict  the 
propagation  characteristics  of  the  beam. 


The  spot  size  is  defined  as  the  radial  distance  from  the  center  to  the 
1/e  point  of  the  Gaussian  beam's  amplitude  distribution.  In  a homogenous 
medium  the  form  of  the  Gaussian  beam  is  [10] 


E(x,y) 


Eo  • WJ  ’ exPHttz-Mz)]  - + 2Rlly)). 


E(x.y)  - Eo 


(62) 


where  r 
R(z) 
w(z) 


k 


n 

x 


(x2  + y2)l/2  = radial  distance, 

radius  of  curvature  of  the  spherical  wave, 

spot  size  at  distance  z, 

minimum  spot  size  (at  z = 0), 

amplitude  coefficient, 

= propagation  constant  of  the  medium, 
refractive  index  of  the  medium, 
optical  wavelength. 
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Figure  3-3.  Foreground  processed  image  of  figure  3-2. 


Figure  3-4.  Convolution  model  of  figure  3-1. 
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Since  we  are  interested  in  intensity  distributions  as  seen  by  the  video 
camera,  the  beam  spot  size,  w,  is  the  1/e2  point  of  the  intensity  distribu- 
tion. The  general  form  of  the  intensity  distribution  [10]  is 

Wr>2  ?r2 

I(x,y)  - E2(x,y)  - E2  . • exp(-^pjjy). 


or 


I(x,y)  « A • exp(-ar2)  (63) 

w 2 

where  A = E2  • 2 i \ and  a = 2/w2(z). 

O W ^ Z ; 

The  A and  a terms  correspond  to  the  parameters  discussed  for  the  Gaussian 
curve-fitting  tool  in  section  2.4. 

Using  the  Gaussian  curve-fit  tool  we  can  measure  the  spot  size,  w,  at  a 
certain  distance.  Subsequent  spot  sizes  can  be  then  be  predicted  by  the  equa- 
tion [10], 

w2(z)  = w 02  • (1  + (64) 


or 


when  ( 


TTW 


X§-)2  » 1. 

o n 


w(z) 


AZ 

ttw  n 
o 


In  the  first  experiment,  as  shown  in  figure  3-5,  we  looked  at  a near- 
field image  of  a beam  that  is  approximately  Gaussian.  We  used  a single-mode 
optical  fiber  that  nominally  has  a 5 pm  core  diameter  and  a numerical  aperture 
of  0.11.  The  fiber  is  illuminated  by  a quartz  halogen  lamp  with  a 851  nm 
narrowband  filter.  The  other  end  of  the  fiber  is  then  magnified  by  successive 
40X  and  10X  microscope  objectives.  The  near-field  image  was  captured  with  a 
vidicon  equipped  with  automatic  gain  control. 


The  digitized  image  of  the  near-field  spot  is  shown  in  figure  3-6  with  a 
resolution  of  256  x 240.  Averaging  four  similar  images  is  necessary  to  reduce 
the  video  noise  resulting  from  the  low-light  conditions.  Figure  3-7  is  the 
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Reticle 


Figure  3-5.  Experimental  setup  for  the  near-field  image  of  an  optical  fiber. 

image  of  a precision  stage  micrometer  placed  in  the  plane  z = 0.  The  lines 
are  2 pm  apart,  so  the  screen  has  a total  horizontal  length  of  24  ym. 

To  obtain  the  beam  spot  radius  we  first  reduced  the  resolution  of  the 
near-field  image  to  64  x 60  using  the  resolution  changing  function.  This 
allowed  for  adjacent  lines  to  be  averaged  to  reduce  possible  errors.  Fig- 
ure 3-8  is  the  intensity  histogram  of  the  near-field  image  and  shows  good  con- 
trast. We  took  the  peak  histogram  intensity  value,  83,  as  the  average  value 
of  the  dark  background  with  0 intensity.  The  intensity  of  the  image  of  the 
near  field  was  then  reduced  by  -83  using  the  biasing  function. 

Figure  3-9  is  the  Gaussian  curve  fit  to  the  near-field  beam.  The  con- 
tinuous line  represents  the  fitted  curve  and  the  dots  represent  the  actual 
data  points.  A noise  cutoff  value  was  n = 20  (sect.  2.4).  The  value  of  n is 
a subjectively  chosen  value  of  the  maximum  noise  shown  in  figure  3-9.  The 
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Figure  3-6.  Near-field  image  of  an  optical  fiber. 


Figure  3-7.  Image  of  scaling  reticle  for  near  field. 
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INI  E N S I T V 

Figure  3-8.  Intensity  histogram  of  near-field  image. 
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Figure  3-9.  Gaussian  curve  fit  of  near-field  image. 
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z =2540pm 


Figure  3-10.  Experimental  setup  for  the  far-field  image  of  an  optical  fiber. 


computed  beam  spot  radius  is  8.2891  pixels  with  64  pixels  across  the  horizon- 
tal screen  length.  The  spot  radius  at  z = 0 is 


,>8.2891  pixels^ 
wo  ^ 64  pixels  ' 


(24  ym)  = 3.108  ym. 


In  the  second  experiment,  shown  in  figure  3-10,  we  looked  at  the  far- 
field  image.  We  used  the  same  single-mode  optical  fiber  and  851  nm  illumina- 
tion. The  far-field  beam  of  the  fiber  at  a distance  z = 2540  ym  was  magnified 
by  a single  lens  and  imaged  onto  the  vidicon  with  automatic  gain  control. 


The  digitized  image  of  the  far-field  beam  spot  is  shown  in  figure  3-11 
with  a resolution  of  256  x 240.  Averaging  eight  similar  images  was  necessary 
to  reduce  the  video  noise  resulting  from  the  low-light  conditions.  Fig- 
ure 3-12  is  the  image  of  a stage  micrometer  placed  in  the  plane  z = 2540  ym. 
Lines  are  10  ym  apart,  so  the  screen  has  a total  horizontal  length  of  about 
950  ym. 


To  obtain  the  beam  spot  radius  we  first  reduced  the  resolution  of  the 
near-field  image  to  64  x 60  using  the  resolution  changing  function.  This 
allows  adjacent  lines  to  be  averaged  to  reduce  possible  errors.  Figure  3-13 
is  the  intensity  histogram  of  the  near-field  image,  and  shows  moderate  con- 
trast. Again,  we  take  the  peak  histogram  intensity  value,  65,  as  the  average 
value  of  the  dark  background  with  0 intensity.  The  intensity  of  the  image  of 
the  near  field  was  therefore  reduced  by  -65  with  the  biasing  function. 
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Figure  3-11.  Far-field  image  of  an  optical  fiber. 


Figure  3-12.  Image  of  scaling  reticle  for  far  field. 
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Figure  3-13.  Intensity  histogram  of  far-field  image. 
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Figure  3-14.  Gaussian  curve  fit  of  far-field  image. 
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Figure  3-14  is  the  Gaussian  curve  fit  to  the  far-field  beam  with  a noise 
cutoff  value  of  32.  The  beam  spot  radius  is  13.629  pixels  with  64  pixels 
across  the  horizontal  screen  length.  The  spot  radius  at  z = 2540  ym  is  there- 
fore 


f 13.629  pixels^  /ncn  \ om  o 
wo  = t~64  -pixel's  • (95°  um)  = 202'3  vm ■ 

We  can  test  the  results  of  the  near-  and  far-field  experiments  by  pre- 
dicting the  far-field  spot  radius  from  the  spot  radius  in  the  near  field. 
Using  the  formula  for  the  propagation  of  the  Gaussian  beam,  eq  (64),  we  see 
that 


,/_N  Xz  (0.851  urn)  (2540  pm)  001  „ 

w(z)  “ -rrwon  (3.1415)  (3.108  pm)(1.0)  221,4  pm’ 

whereas  we  have  measured  a value  of  202.3  pm.  The  relative  error  is  thus 


% error  = - . 10q%  ~ 9%. 

202.3  pm 


The  near-  and  far-field  images  produce  consistent  results. 

The  error  incurred  in  this  set  of  preliminary  experiments  is  a result  of 
many  contributions.  The  small  distance  z was  hard  to  measure  with  great 
accuracy  but  can  be  increased  through  the  use  of  longer  focal  length  lenses. 
The  recti  1 inearity  and,  especially,  the  linearity  of  responsivity  of  the  vidi- 
con  camera  are  questionable.  These  could  be  improved  greatly  with  a CCD 
(charge  coupled  device)  array  camera  but  light  from  the  beam  signal  was  too 
low  for  our  CCD  camera.  This  may  be  improved  with  a broader  band  interference 
filter  or  a laser  source.  Distortion  in  the  lenses  is  also  subject  to  scru- 
tiny. 


4.  Conclusion 

4.1  Conclusions  and  Future  Applications 

The  purpose  of  this  thesis  was  to  develop  and  test  image  processing  soft- 
ware designed  for  optical  engineering  applications.  Image  processing  is  a 
relatively  new  and  rapidly  advancing  field.  Present  applications  of  image 
processing  include  aerial  photography  analysis,  medical  analysis,  photographic 
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processing,  robotic  image  recognition,  and  inspection,  to  name  a few.  Most  of 
the  software  and  systems  available  today  are  designed  with  universal  applica- 
tions in  mind  to  meet  the  widest  sector  of  potential  customers  in  this  new 
area.  The  software  in  this  thesis  is  partly  universal  in  nature,  but  has  been 
written  with  optical  engineering  applications  in  mind. 

Image  processing  has  many  strengths  which  are  advantageous  to  optical 
engineering.  By  using  software  and  mathematically  modeling  optical  devices  we 
can  simulate  optical  processing  without  the  actual  devices.  This  is  important 
in  the  design  and  simulation  of  electrical  and  electronic  components  in  very 
large  scale  integrated  (VLSI)  electronic  chips.  Simulating  optical  processing 
may  prove  to  be  even  more  important. 

Since  image  processing  deals  with  digital  images  we  can  process  such 
images  with  great  numerical  accuracy.  Real  optical  components  are  also  image 
processors  but  are  quite  costly  to  make  with  high  quality  and  are  time  consum- 
ing to  implement.  Image  processing  offers  a cheaper  and  faster  alternative 
with  greater  possibilities.  The  image  processing  functions  described  in 
section  2 process  images  without  coherent  optics.  With  more  sophisticated 
hardware  and  algorithms  the  software  may  actually  do  a better  job  than  real 
optical  devices.  Computer  processing  can  also  perform  functions  that  cannot 
be  performed  optically,  and  sometimes  these  can  be  used  to  advantage  in  coher- 
ent optical  processing. 

As  shown  in  section  3.2,  we  can  also  use  image  processing  in  fiber  and 
integrated  optics;  other  optical  devices  can  also  be  analyzed  with  computer 
processing.  This  is  important  in  both  the  characteri zation  and  testing  of 
experimental  and  newly  manufactured  optical  devices.  The  output  of  whole 
optical  systems  can  also  be  analyzed. 

The  software  in  this  thesis  serves  as  an  introduction  to  image  processing 
for  optical  engineering.  Possible  future  directions  for  a project  of  this 
type  involve  a system  with  greater  resolution  and  computing  capacity.  Today's 
digital  design  is  heavily  involved  in  working  to  create  faster  and  more  spe- 
cialized hardware  to  meet  these  goals.  With  greater  computing  power  more  com- 
plex and  challenging  problems  can  be  investigated.  There  is  also  the  need  for 
better  input  and  output  devices  for  image  processing  systems.  High  resolution 
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digitizing  cameras  and  scanning  devices  are  very  expensive  and  limited.  Simi- 
larly, high  resolution  displays  and  printers  are  advancing  but  still  have  a 
long  way  to  go.  In  particular,  a cheap,  fast,  high  resolution  coherent  dis- 
play would  be  most  useful.  Liquid  crystal  televisions  are  cheap  and  avail- 
able, but  have  low  resolution,  contrast,  and  speed.  Ferroelectric  crystal 
displays  are  more  promising  in  these  areas  but  are  expensive  and  not  fully 
developed  at  present.  For  both  input  and  output  devices  there  is  a need  for 
recti  1 inearity  and  linearity  of  responsivity  to  ensure  high  image  fidelity  and 
measurement  accuracy. 
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Appendix  A.  User's  Manual 


1.  Introduction 

1.1  Software  Description 

IMPROC  is  an  image  processing  software  designed  primarily  for  optical 
engineering  applications.  It  contains  a unique  set  of  image  processing  tools 
not  found  in  other  current  software. 

IMPROC  is  generic  in  the  sense  that  it  works  on  ordinary  MS-DOS  data 
files.  IMPROC  is  designed  to  be  used  in  conjunction  with  a video  image  digit- 
izer or  frame-grabber.  Picture  files  created  by  the  frame-grabber  are  used  as 
the  input  to  IMPROC.  IMPROC  can  manipulate  and  display  information  about 
these  files  using  the  various  image  processing  tools.  These  altered  files  can 
be  stored  again  as  picture  files,  which  can  be  subsequently  displayed  by  the 
frame  grabber. 

2.  Program  Information 

2.1  Picture  File  Specifications 

Picture  files  are  the  main  means  of  storage  of  a digitized  image.  Each 
pixel  of  the  digitized  image  is  represented,  by  one  byte.  This  allows  for  the 
capability  of  up  to  256  levels  of  intensity  (gray  scale)  per  pixel.  The  pic- 
ture file  consists  of  a linear  array  of  pixels  which  represent  the  two-dimen- 
sional image.  The  file  has  a resolution  format  which  describes  the  vertical 
and  horizontal  resolution  of  the  image.  The  valid  resolution  formats  allowed 
by  IMPROC  are 


Resolution  format 
(vert,  by  horiz.) 

File  s i ze 
(bytes) 

16x15 

240 

32x30 

960 

64x60 

3,840 

128x120 

15,360  (15K) 

256x240 

61,440  (60K) 
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I MPROC  allows  the  user  to  choose  between  detail  and  computational  speed  by 
changing  resolution  formats.  Input  picture  files  with  resolution  formats  not 
listed  above  are  illegal.  Illegal  files  cannot  be  used  and  will  prompt  an 
error  when  attempted  to  be  read  by  IMPROC. 

Also  throughout  the  software  certain  conventions  about  directional 
indices  and  implied  center  points  are  assumed. 


Res.  Format 

Horiz.  Indices 

Vert.  Indices 

Center 

16x15 

0 , 1 , . . . , 15 

0,1,. ..,14 

(8,7) 

32x30 

0,1,. ..,31 

0,1,. ..,29 

(15,14) 

64x60 

0,1 , . . . ,63 

0, 1, ... ,59 

(31,29) 

128x120 

0,1,. ..,127 

0,1,. ..,119 

(63,59) 

256x240 

0,1, .. . ,255 

0,1,. ..,239 

(127,119) 

These  conventions  are  necessary  for  correct  interpretation  and  manipulation  in 
center  dependent  functions  such  as  Fourier  transform  and  magnif ication. 

2.2  Program  Execution 

The  IMPROC  program  is  initiated  by  simply  typing  IMPROC.  This  calls  the 
MS-DOS  batch  file  IMPROC.BAT,  which  is  the  heart  of  the  IMPROC  software.  This 
file  controls  the  program  flow  of  the  three  IMPROC  execution  files, 

IMPR0C#1 .EXE,  IMPR0C#2.BAS,  and  IMPR0C#3.EXE. 

IMPR0C#1.EXE  is  an  executable  compiled  C program  which  contains  all  the 
picture  file  manipulation  code  and  is  by  far  the  largest  of  the  three.  This 
portion  of  the  code  is  written  in  C to  take  advantage  of  its  speed  and  program 
flow,  which  are  necessary  for  the  long  and  complex  image  processing  algo- 
rithms. 

IMPR0C#2.BAS  is  an  interpreted  BASIC  program  used  for  two  dimensional 
graphing  on  the  graphics  monitor.  IMPR0C#3.EXE  is  an  executable  compiled 
BASIC  program  used  for  three-dimensional  plotting  on  the  graphics  monitor. 
These  two  programs  were  written  in  BASIC  because  of  the  lack  of  convenient 
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graphic  capabilities  in  the  C language.  The  second  BASIC  file,  IMPR0C#3.EXE, 
is  compiled  to  speed  up  the  tedious  process  of  plotting. 

Information  is  passed  among  the  three  programs  by  temporary  data  files, 
IMPR0C#1 .PIC,  IMPR0C#2 .DAT,  and  IMPR0C#3.DAT.  IMPR0C#1 .PIC  is  the  current  pic- 
ture being  used  by  IMPR0C#1.EXE  and  is  stored  when  a call  to  IMPR0C#2.BAS  or 
IMPR0C#3.EXE  is  made.  IMPR0C#2.DAT  is  the  graphing  data  created  by 
I MPR0C#1 .EXE  and  is  used  by  the  graphing  program  IMPR0C#2.BAS.  IMPR0C#3.DAT 
is  the  plotting  data  created  by  IMPR0C#1.EXE  and  is  used  by  the  plotting  pro- 
gram IMPR0C#3.EXE. 

2.3  Program  Compilation 

To  make  changes,  additions,  or  corrections  in  the  software,  we  need  to 
know  how  to  compile  the  program  source  code.  The  program  is  written  C using 
Lattice  C Rev.  2.14  [11].  The  program  should  be  compiled  under  the  "d"  model 
which  stresses  a moderately  sized  program  which  works  with  a large  amount  of 
memory.  The  batch  file,  CCDFFT.BAT,  controls  the  compilation  process 

lei  %1  -mD 
1 c2  %1 

link  cd+fft87+%l ,%1 ,nul ,lcmd+lcd 

To  compile  we  simply  need  to  type  CCDFFT  IMPR0C#1  and  return.  The  machine 
code  Fourier  transform  program,  FFT87.C  (Rapid  Imaging  Software  [11]),  is 
automatically  linked.  Compilation  takes  the  source  file,  IMPR0C#1.C,  and 
generates  the  executable  file  IMPR0C#1 .EXE. 

3.  Commands 

3.1  Image  Commands 

3.1.1  Average 

COMMAND:  i [RETURN]  a [RETURN] 

FUNCTION:  p(x)  = ( p(x)  + input(x)  ) / 2 

DESCRIPTION:  Averages  the  individual  pixels  of  the  current  picture  in 

memory  with  an  input  picture  file.  The  result  becomes  the  current  picture. 
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3.1.2  Bias 


COMMAND:  i [Return]  b [Return] 

FUNCTION:  p(x)  = p(x)  + bias 

DESCRIPTION:  Adds  an  input  integer  value  to  each  pixel.  Valid  range  is 

-255  to  255.  Biased  values  >255  are  set  to  255  and  values  <0  are  set  to  0. 

OPTION(S):  The  auto-bias  option  computes  a histogram  and  automatically 

biases  (negatively)  with  the  peak  relative  frequency  pixel  value. 

3.1.3  Convolution 

COMMAND:  i [RETURN]  c [RETURN] 

FUNCTION:  p(x)  = cross  correlation  with  input  matrix 

DESCRIPTION:  Computes  the  cross  correlation  of  the  neighbors  of  each 

pixel  with  an  input  matrix  and  sets  the  value  to  that  pixel. 

The  input  matrix  has  two  dimensions,  each  of  which  must  have  an  odd 
integer  length.  The  matrix  is  centered  on  the  current  pixel  and  the  elements 
of  the  matrix  must  be  integers.  A matrix  normalization  value  is  automatically 
calculated  to  set  the  area  of  the  matrix  elements  to  1.  If  the  matrix  has 
zero  area  the  normalization  constant  is  set  to  1. 

During  convolution  the  elements  of  the  matrix  are  multiplied  with  their 
respective  pixel  neighbors.  The  sum  of  these  products  is  then  normalized  and 
biased  (optional).  Pixels  off  the  edge  will  cause  erroneous  convolutions  and 
should  be  ignored.  Pixel  values  created  which  are  >255  are  set  to  255  and 
values  <0  are  set  to  0. 

OPTION(S):  The  normalization  value  may  be  set  manually  but  not  to  zero. 

Also  a bias  value  may  be  added  after  the  convolution  of  each  pixel. 
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3.1.4  Demagnify 


COMMAND:  i [RETURN]  d [RETURN] 

FUNCTION:  Demagnify  picture  by  2X 

DESCRIPTION:  Shrinks  the  current  picture  by  a factor  of  2 in  each  direc- 

tion and  places  it  in  the  center.  This  process  results  in  a 4X  loss  of  infor- 
mation as  averaging  of  four  pixels  is  used  to  create  each  demagnified  pixel. 
The  pixels  of  the  blank  area  around  the  demagnified  picture  are  set  to  0. 

3.1.5  Enhance  Contrast 

COMMAND:  i [RETURN]  e [RETURN] 

FUNCTION:  p(x)  = (255/ (hi ghest-1 owest ) ) * (p(x)-lowest) 

DESCRIPTION:  Increases  the  contrast  over  a specified  range  [lowest,  high- 

est] to  full  contrast  [0,255].  Pixels  values  >255  are  set  to  255  and  values 
<0  are  set  to  0. 

OPTION(S):  The  autoscale  option  automatically  finds  the  lowest  and  high- 

est values  to  enhance  the  contrast  to  full  scale  without  clipping. 

3.1.6  Fourier  Transform 

COMMAND:  i [RETURN]  f [RETURN] 

FUNCTION:  creates  a 2-D  FFT  of  picture. 

DESCRIPTION:  Creates  a two-dimensional  Fourier  transform  using  a one- 

dimensional fast  Fourier  transform  algorithm  (FFT).  The  transform  assumes 
that  the  picture  is  an  intensity  distribution  made  up  of  real  components.  The 
square  root  of  each  pixel  is  taken  to  represent  the  light  amplitude.  First 
the  function  processes  the  picture  in  the  vertical  direction  to  obtain  the 
vertical  transforms.  Second  the  function  processes  in  the  horizontal  direc- 
tion to  obtain  the  horizontal  transforms.  The  result  is  taken  to  be  the 
square  of  the  real  plus  the  square  of  the  imaginary  components  which  represent 
the  Fourier  transform  intensity  distribution. 

OPTION(S):  The  user  can  make  the  choice  of  obtaining  either  the  forward 

or  inverse  Fourier  transform. 

The  user  can  also  select  an  appropriate  power-of-2  Fourier  plane  magnifi- 
cation to  get  a reasonably  sized  output  transform  relative  to  the  input 
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picture.  Fourier  plane  magnification  results  in  the  use  of  larger  FFT  arrays 
and  subsequently  requires  more  computation. 

Finally  the  Fourier  transform  must  be  scaled  to  fit  with  the  pixel  limits 
[0,255].  The  maximum  point  of  the  transform  is  taken  to  be  255.  Using  a 
scaling  factor  provides  that  the  output  can  be  scaled  to  observe  finer  detail 
not  found  in  the  dynamic  range  of  the  pixels.  The  height  of  the  intensity 
pattern  can  then  be  taken  to  be  the  maximum  point  divided  by  the  scale  factor. 

3.1.7  Magnify 

COMMAND:  i [RETURN]  m [RETURN] 

FUNCTION:  Magnify  picture  by  2X. 

DESCRIPTION:  Enlarges  the  center  of  the  current  picture  by  a factor  of  2 

in  each  direction.  This  process  results  in  a 4X  loss  of  information  as  the 
pixels  outside  the  center  to  be  magnified  are  lost.  Center  pixels  are  repli- 
cated four  times  each  to  create  the  new  picture. 

3.1.8  Negative 

COMMAND:  i [RETURN]  n [RETURN] 

FUNCTION:  p(x)  = 255  - p(x) 

DESCRIPTION:  Creates  a reversed  contrast  image. 

3.1.9  Pratt's  Noise  Cleaning  Algorithm 

COMMAND:  i [RETURN]  p [RETURN] 

FUNCTION:  Reduces  noise  from  picture. 

DESCRIPTION:  Cleans  noise  generated  by  a video  camera  or  image  process- 

ing. The  noise  threshold  is  set  by  the  input  value  of  e.  If  the  absolute 
value  of  the  difference  between  a pixel  value  and  the  average  of  its  eight 
neighbors  is  greater  than  e,  then  the  pixel  value  is  set  to  the  average  of  its 
nei ghbors . 
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3.1.10  Resolution  Change 


COMMAND:  i [RETURN]  r [RETURN] 

FUNCTION:  Changes  resolution  format. 

DESCRIPTION:  Reduces  or  expands  the  vertical  by  horizontal  resolution 

format  of  the  current  picture.  In  resolution  reduction  the  pixels  are  aver- 
aged; this  results  in  an  information  loss,  and  the  new  picture  has  fewer 
pixels.  In  resolution  expansion  the  pixels  are  replicated  with  no  information 
gain  although  the  new  picture  has  more  pixels.  If  the  resolution  is  chosen  to 
remain  the  same  the  picture  is  unchanged. 

3.1.11  Threshold 

COMMAND:  i [RETURN]  t [RETURN] 

FUNCTION:  if  ( p(x)  >=  lower  bound  ) and  ( p(x)  <=  upper  bound  ) then 

p(x)  = new  value 

DESCRIPTION:  Sets  each  pixel  to  a new  value  if  that  pixel  is  within  the 

desired  range.  The  range  is  determined  by  the  upper  and  lower  bound. 

3.2  Display  Commands 
3.2.1  Gaussian  Curve  Fit 

COMMAND:  d [RETURN]  g [RETURN] 

FUNCTION:  Fits  a Gaussian  curve  to  a picture. 

DESCRIPTION:  Calculates  the  centroid  or  first  moment  of  a Gaussian  image. 

A horizontal  or  vertical  line  that  passes  through  the  centroid  may  be 
selected.  A Gaussian  curve  fit  is  performed  using  a least  squares  fit  to  the 
linearized  data  of  the  line. 

The  1/e2,  5%,  4%,  and  3%  points  are  calculated  from  the  parameters  of  the 
fitted  Gaussian  curve.  Also  a plot  of  the  selected  line  is  generated  with  the 
Gaussian  curve  shown  as  a continuous  line  and  the  actual  data  points  repre- 
sented by  dots. 

OPTION(S):  We  must  enter  a minimum  noise  cutoff  value.  Pixel  values  less 

than  or  equal  to  this  value  are  ignored  in  calculating  the  centroid  and 
Gaussian  curve  as  these  values  are  taken  to  be  noise.  We  can  also  select 
either  a horizontal  or  vertical  fit. 
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3.2.2  Multimode  Curve  Fit  (g-Profile) 


COMMAND:  d [RETURN]  m [RETURN] 

DESCRIPTION:  Calculates  the  centroid  of  the  multimode  profile.  A hori- 

zontal or  vertical  line  that  passes  through  the  centroid  may  be  selected.  A 
g-Profile  fit  is  performed  using  a least-squares  fit  on  the  data  in  the  10 
percent  to  80  percent  range  (EIA  standard)  of  Imax  which  is  the  maximum  pixel 
value  of  the  digitized  image.  The  parameters  a (core  radius)  and  g are 
cal cul ated. 

OPTION(S):  We  must  enter  a minimum  noise  cutoff  value.  Pixel  values  less 

than  that  or  equal  to  this  value  are  ignored  in  the  calculation  of  the  cen- 
troid. We  can  also  enter  a radius  of  exclusion  which  omits  all  data  points 
within  a certain  radius  from  the  center. 

A horizontal  or  vertical  fit  may  be  selected. 

3.2.3  Intensity  Histogram 

v 

COMMAND:  d [RETURN]  h [RETURN] 

DESCRIPTION:  Calculates  the  relative  frequency  versus  pixel  value  over 

the  entire  picture  and  creates  a graph. 

OPTION (S):  The  output  graph  can  be  displayed  in  either  the  line  or  bar 

graph  mode.  If  the  bar  graph  mode  is  chosen  the  bar  width  may  be  specified. 

A one-line  message  may  be  entered  to  be  displayed  on  the  output  graph. 

3.2.4  Line  Intensity  Scan  (2-D) 

COMMAND:  d [RETURN]  2 [RETURN] 

DESCRIPTION:  Displays  intensity  as  a function  of  position  of  a line 

across  the  picture. 

OPTIONS:  The  user  can  use  either  a horizontal  or  vertical  scan  across  the 

picture.  If  a horizontal  scan  is  selected  the  user  must  input  the  vertical 
index  or  position  of  the  scan.  If  a vertical  scan  is  selected  the  user  must 
input  the  horizontal  index  or  position  of  the  scan.  The  output  graph  can  be 
displayed  in  either  the  line  or  bar  graph  mode. 

A one-line  message  may  be  entered  to  be  displayed  on  the  output  graph. 
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3.2.5  3-D  Intensity  Plot 


COMMAND:  d [RETURN]  3 [RETURN] 

DESCRIPTION:  Creates  a three-dimensional  intensity  plot  of  the  picture. 

The  plot  can  have  a maximum  of  60  horizontal  lines  and  256  points  per  line. 
Since  the  function  is  limited  to  only  60  horizontal  lines  all  of  the  informa- 
tion of  the  highest  resolution  pictures  (128x120  or  256x240)  cannot  be  dis- 
played. 

OPTION(S):  In  the  highest  resolution  pictures  (128x120  or  256x240)  the 

user  may  select  a 4X  magnification  around  the  center  point. 

Autoscaling  can  be  used  to  find  the  lowest  and  highest  pixels  contained 
in  the  picture  to  ensure  that  the  plot  has  the  appropriate  height. 

We  can  chose  between  a single-line  or  a cross-hatch  type  of  3-D  graph. 
Also  a one-line  message  may  be  entered  to  be  displayed  on  the  output  plot. 

3.3  File  Commands 

3.3.1  Read 

COMMAND:  r [RETURN] 

DESCRIPTION:  Reads  a picture  data  file  into  memory.  To  read  the  file  the 

user  must  input  the  complete  picture  file  name  (and  path  if  necessary).  After 
the  file  is  read  the  number  of  bytes  read  and  resolution  format  of  the  picture 
are  displayed.  On  an  illegal  resolution  format,  an  error  will  occur  (see  pic- 
ture file  specifications,  section  2.2). 

3.3.2  Write 

COMMAND:  w [RETURN] 

DESCRIPTION:  Writes  a picture  in  memory  to  a file.  To  write  the  file  the 

user  must  input  the  complete  picture  file  name  (and  path  if  necessary). 
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Appendix  B.  Source  Code 


The  source  code  for  this  software  is  made  up  of  four  executable  files. 

The  first  file,  IMPROC.BAT  is  the  main  batch  file  responsible  for  controlling 
the  other  three  files.  The  second  file,  IMPR0C#1 . EXE , is  an  executable  com- 
piled C file  of  the  image  processing  program  IMPR0C#1.C.  The  third  file, 
IMPR0C#2.BAS,  is  a BASIC  program  which  creates  the  two-dimensional  graphs. 

And  the  fourth  file,  IMPR0C#3.EXE,  is  an  executable  compiled  BASIC  file  of  the 
program  IMPR0C#3.BAS  which  creates  the  three-dimensional  plots. 

There  are  also  other  necessary  files  which  are  not  listed  in  this  appen- 
dix. The  file,  INIT_FFT.COM  (Rapid  Image  Software,  Tijeras,  N.M.),  is  a 
machine  code  FFT  program  and  must  be  executed  before  running  the  software. 

The  file,  GRAPHICS.COM  (IBM  Corp.)  must  also  be  executed  before  running  the 
software  for  proper  printing  of  graphics.  The  files,  BASICA.COM  and 
BASRUN.EXE  must  be  present  for  the  execution  of  IMPR0C#2.BAS  and  IMPR0C#3.EXE 
[11]. 

IMPROC . BAT 


echo  off 

if  exist  improc#l . pic  erase  improc#l.pic 
if  exist  improc#2 . dat  erase  improc#2 . dat 
if  exist  improc#3 . dat  erase  improc#3 . dat 
: loop 
mode  mono 
els 

improc#l  =32768 

if  exist  improc#2.dat  mode  co80 
if  exist  improc#2 . dat  basica  improc#2 
if  exist  improc#3.dat  mode  co80 
if  exist  improc#3.dat  improc#3 
if  not  exist  improc#l.pic  goto  end 
goto  loop 
: end 

mode  mono 
els 
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IMPROC#! .C 


#include  <stdio.h>  /* 
#include  <fcntl.h>  /* 
# include  <math.h>  /* 
#include  <limits.h>  /* 

#def ine  MAXLINE  40  /* 
#def ine  MAXMAT  50  /* 
#def ine  PSIZESIZE  4 /* 
#def ine  MAXPSIZE  61440  /* 
#def ine  FBYTES  240  /* 
#def ine  MAXFFT  4096  /* 
#def ine  EIAMMMAX  0.8  /* 
#def ine  EIAMMMIN  0.1  /* 


standard  input/output  routines  */ 
memory /hardware  functions  */ 
math  routines  */ 
standard  limits  */ 

# of  char,  read  as  screen  input  */ 
#-l  of  non-zero  convolution  elements 

# of  elements  in  psize  array  */ 
maximum  picture  size  (60K)  */ 
file  10  block  size  */ 
maximum  FFT  array  length  */ 

EIA  multi-mode  maximum  point  */ 

EIA  multi -mode  minimum  point  */ 


*/ 


main( ) 

{ 

char  line [MAXLINE] , *malloc (), *p ; 

/*  line:  array  for  user  input  */ 

/*  malloc():  memory  allocation  function  */ 

/*  p[]:  single  dimension  picture  array  */ 
unsigned  psize [ PSIZESIZE ] ; 

/*  psize[]:  picture  size  attributes  array 
psize[0]:  total  picture  array  length 
psizefl]:  vertical  picture  length 
psize[2]:  horizontal  picture  length 

psize[3]:  picture  resolution  format  number  (1..5)  */ 
int  *exit , quit , ret ; 

/*  exit:  program  exit  flag  */ 


p=malloc (MAXPSIZE) ; /*  picture  array  memory  allocation  */ 

if  (p==NULL)  printf( "*****  malloc  error  *****\n"); 
init_pic(p , psize) ; 

/*  reset  interupted  picture  array  from  a previous  program  to  execute  a 
graphics  operation  (BASIC)  */ 


quit=0 ; 

while  (quit==0)  { 
printf ("\n") ; 
printf("Main  Menu:\n"); 

printf("\td  = display  functions  (graphics)\n" ) ; 
printf("\ti  = image  processing  functions\n" ) ; 
printf("\tr  = read  a file\n"); 
printf ("\tw  = write  a file\n"); 
printf("\tq  = quit\n"); 


printf ( "\n===>  ");  ret=getline ( line) ; 
if  (ret!=0)  { 

switch  (line[0])  { 
case ' D ' : 

case'd':  display (p , psize , exit) ; 

if  (*exit==l)  { save_pic (p , psize) ; quit=l ; } ; 

/*  save  current  picture  before  exiting  program  */ 
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break; 


case ' I ' 
case ' i ' 
case ' R' 
case ' r ' 
case ' W' 
case ' w' 
case'Q' 
case ' q ' 
default 


image (&p , psize) ; break; 
read_pic (p , psize)  ; break; 
write_pic (p , psize) ; break; 
quit=l ; break; 

printf ( " %c*****  invalid  command 


*****\n" , 7 ) ; break ; 


) ; 

ret=free(p);  /*  free  picture  memory  */ 


display (p , psize , exit) 
char  *p; 

unsigned  *psize; 
int  *exit; 

{ 

char  line [MAXLINE] ; 
int  ret; 


printf ("\n") ; 

printf ( "Display  Menu:\n"); 
printf("\tg  = gaussian  curve  fit\n"); 
printf("\th  = histogram  of  intensity  graph\n"); 
printf("\tm  = multi-mode  curve  fit  (g-prof ile)\n" ) ; 
printf("\t2  = 2-D  line  intensity  graph\n"); 
printf("\t3  = 3-D  intensity  plot\n"); 
printf("\tq  = quit\n"); 


=>  ");  ret=getline ( line) ; 


printf  ( "\n= 

*exit=0 ; 
if  (ret!=0)  { 

switch  (line[0])  { 

/*  exit  will  end  the  program  to  execute  graphics  program  (BASIC)  */ 
case ' G ' 

gaussian(p , psize) ; *exit=l ; break; 


case ' g' 
case ' H' 
case 'h' 
case 'M' 
case ' m' 
case ' 2 ' 
case ' 3 ' 
case ' Q' 
case ' q ' 
default 


histogram(p , psize)  ; *exit=l ; break; 

multi_mode (p , psize) ; *exit=l ; break; 
graph_2d(p , psize) ; *exit=l ; break; 
plot_3d(p , psize) ; *exit=l ; break; 

break ; 

printf ("%c*****  invalid  command  *****\n" , 7) ; break; 


) ; 
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image (pp , psize) 
char  **pp ; 
unsigned  *psize; 

{ 

char  i ine[ MAXLINE ] ,*p; 
int  mv[MAXMAT] , ret ; 

/*  mv[]:  convolution  matrix  value  array  */ 
unsigned  ma[MAXMAT] ; 

/*  ma[]:  convolution  matrix  relative  address  array  */ 


printf ("\n") ; 

printf ("Image  Menu:\n"); 

printf("\ta  = average  with  another  file\n"); 

printf("\tb  = bias\n"); 

printf("\tc  = convolve\n" ) ; 

printf("\td  = demagnif ication\n" ) ; 

printf("\te  = enhance  contrast\n" ) ; 

printf("\tf  = fourier  transform  (2-D  FFT)\n"); 

printf("\tm  = magnif ication\n" ) ; 

printf("\tn  = negative  (reverse  contrast) \n" ) ; 

printf("\tp  = Pratt's  noise  cleaning  algorithm\n" ) ; 

printf ("\tr  = resolution  change\n"); 

printf("\tt  = threshold\n" ) ; 

printf("\tq  = quit\n"); 

printf ( "\n======>  ");  ret=getline(line) ; 

if  (ret !=0)  { 

P=*PP ; 

switch  (line[0])  { 
case 'A' 

average(p , psize) ; break; 


case ' a' 
case ' B ' 
case 'b ' 
case ' C ' 
case ' c ' 
case ' D ' 
case ' d' 
case ' E' 
case ' e ' 
case ' F' 
case ' f ' 
case ' M' 
case ' m' 
case ' N' 
case ' n' 
case ' P ' 
case ' p ' 
case ' R' 
case ' r ' 
case ' T ' 
case ' t ' 
case ' Q ' 
case ' q ' 
default 


bias (p , psize) ; break; 

convolve(pp, psize, ma,mv, 1,0);  break; 

demag(pp , psize) ; break; 

enhance (p , psize) ; break; 

fourier (p , psize) ; break; 

mag(pp , psize) ; break; 

negative(p , psize) ; break; 

pratt(pp , psize ,ma,mv) ; break; 

resolution(pp , psize) ; break; 

threshold(p , psize) ; break; 

break; 

printf (" %c*****  invalid  command  *****\n" , 7 


) ; break 
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) 


) ; 


); 


average (p , psize) 
char  *p ; 

unsigned  *psize; 

( 

char  *malloc ( ) , *temp ; 
int  ret; 

unsigned  psizeO , tempsize [PSIZESIZE] ,u; 
temp=malloc (MAXPSIZE) ; 

if  (temp==NULL)  printf ( "%c*****  malloc  error  *****\n",7); 
else  { 

read_pic ( temp , tempsize) ; 
if  ( tempsize [ 0 ]! =psize [ 0 ] ) 

printf (" %cResolution  of  input  does  not  match  current  file.", 7); 
else  { 

printf ( "Averaging  with  input  file. . . \n") ; 
psizeO=psize [ 0 ] ; 

for  (u=0;  u<psizeO;  ++u)  p [u]=(p [u]+temp [u]+l)/2 ; 

); 

ret=f ree ( temp) ; 

) ; 

} 

bias (p , psize) 
char  *p ; 

unsigned  Upsize; 

{ 

int  biasvalue , exit , mode , peak , value ; 
unsigned  h256 [ 256 ], psizeO , u; 

printf ( "Select  mode  (0=manual , l=auto-bias)  [0]:  ");  mode=getnum( ) ; 
if  (mode!=l)  mode=0 ; 
psizeO=psize [ 0 ] ; 

if  (mode==0)  { 
exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  bias  value  [ - 255 , . . . , 255 ] : ");  biasvalue=getnum( ) 
if  ( (biasvalue>=- 255 ) &&(biasvalue<=255 ) ) exit=l ; 

}; 

); 

if  (mode==l)  { 

printf ("Computing  Histogram. . . \n") ; 
for  (u=0 ; u<256 ; ++u)  h256[u]=0; 
for  (u=0;  u<psize0;  ++u)  ++h256 [p [u] ] ; 
peak=0 ; 

for  (u=0;  u<256;  ++u)  { 
if  (h256 [u] >peak)  { 

peak=h256 [u] ; biasvalue=u; 
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) ; 

); 

printf("Peak  value  is  %d\n" ,biasvalue) ; biasvalue=-biasvalue ; 

) ; 

printf ( "Adding  bias  of  %d. . . \n" ,biasvalue) ; 
for  (u=0;  u<psizeO;  ++u)  { 
value=p [u]+biasvalue ; 
if  (value<0)  p[u]=0; 
else  { 

if  (value>255)  p[u]=255; 
else  p[u]=value; 

} ; 

) ; 

} 

calc_res (psize) 
unsigned  Upsize; 

{ 

switch  ( (int) (psize [0]/2) ) { 


case 

120: 

psize 

1]=  15 

psize 

2]=  16 

psize 

[ 3 ] =1 

break; 

case 

480: 

psize 

1]=  30 

psize 

2]=  32 

psize 

[ 3 ]=2 

break ; 

case 

1920: 

psize 

1]=  60 

psize 

2 ] = 64 

psize 

[ 3 ] =3 

break; 

case 

7680: 

psize 

1 ] =120 

psize 

2 ] =128 

psize 

[ 3 ] =4 

break ; 

case 

30720: 

psize 

1 ] =240 

psize 

2 ] =256 

psize 

[ 3 ]=5 

break ; 

default 

printf (" %c*****  invalid  file  size 
psize[0]=0;  psize[l]=0;  psize[2]= 

■k'k'k'k't 

0 ; ps: 

An",  7); 
Lze [ 3 ] =0 

break ; 

}; 

} 

convolve (pp , psize , ma , mv , option , epsilon) 
char  **pp ; 

unsigned  *psize , *ma ; 
int  *mv, option, epsilon; 

/*  option  1 is  a user  inputed  convolution,  option  2 is  for  Pratt's  algorithm  */ 
/*  epsilon:  threshold  value  from  Pratt's  algorithm  */ 

{ 

char  *malloc ( ) , *p , *temp ; 
int  bias , diff, i ,n,norm, sum; 
unsigned  psizeO ,psize2 ,u,v; 

p=*pp ; temp=malloc(MAXPSIZE) ; 

if  (temp==NULL)  printf ( "*****  malloc  error  *****\n"); 
else  { 

psizeO=psize [0] ; psize2=psize [ 2 ] ; 
bias=0 ; 

if  (option==l)  { 

input_matrix(psize , ma , mv) ; 
printf ("\n") ; 

printf ( "Enter  bias  value  [0]:  ");  bias=getnum( ) ; 

}; 

n=ma[0];  norm=0 ; 
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for  (i-1;  i<=n;  ++i)  norm+=mv[ i] ; 
if  (norm==0)  norm=l ; 
if  (option==l)  { 

printf ("Enter  normalization  value  [%d]:  ",norm);  i=getnum(); 
if  (i!=0)  norm=i; 

) ; 

printf ("Bias  = %d  Normalization  = %d\n" , bias , norm) ; 

printf("  0%%  completed\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ) ; 
u=0 ; 

while  (u<psizeO)  { 

for  (v=0;  v<psize2;  ++v)  { 
sum=0 ; 

for  (i=l;  i<=n;  ++i) 

if  (u+ma[i]<psizeO)  sum+=mv[ i ] *( int)p [u+ma [ i ] ] ; 
sum=sum/norm ; 
if  (option==l)  { 
sum+=bias ; 

if  (sum<0)  temp[u]=0; 

else  if  (sum>255)  temp[u]=255; 

else  temp[u]=sum; 

} 

else  { 

dif f=sum- ( int)p [u] ; 
if  (diff<0)  diff=( -diff ) ; 
if  (diff<epsilon)  temp [u]=p[u] ; 
else  temp[u]=sum; 

); 

++u  ; 

}; 

printf ( " %3 . Of\b\b\b" , ( 100 . 0*u) /psizeO) ; 

}; 

printf ( "\n" ) ; 

*pp=temp ; i=free(p); 

) ; 

) 

demag (pp ,psize) 
char  **pp ; 
unsigned  Upsize; 

{ 

char  *malloc ( ) , *p , *temp ; 
int  ret, sum; 

unsigned  psizeO, psizel ,psize2,px,py, tempx , tempxs tart , tempxf inish , 
tempy , tempys tart , tempyf inish , u ; 

printf  ( "Demagnifying  by  2X.  . .\n") 
p=*pp;  temp=malloc (MAXPSIZE) ; 

if  (temp==NULL)  printf  (" %c**'**'*  malloc  error  ■****"*\n"  , 7 ) ; 
else  { 

psizeO=psize [0] ; psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
for  (u=0;  u<psize0;  ++u)  temp[u]=0; 

tempxs tar t=psize2/4 ; tempxf inish=( 3*p size 2) /4 ; 
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tempystart=(psizel/4)*psize2 ; tempyf inish=( (3*psizel ) /4)*psize2 ; 
py=0; 

for  ( tempy=tempystart ; tempy<tempyf inish ; tempy+=psize2)  ( 
px=0; 

for  ( tempx=tempxstart ; tempxCtempxf inish ; ++tempx)  { 

sum=p [px+py ] +p [px+l+py ] +p [px+py+psize2 ] +p [ px+l+py+psize2 ] ; 
temp [ tempx+tempy ] = ( sum+2 ) /4 ; 
px+=2 ; 

} ; 

py+=(2*psize2) ; 


*pp=temp ; ret=free(p); 

); 

} 

enhance (p , psize) 
char  *p ; 

unsigned  Upsize; 

{ 

int  exit , low, high, mode ; 
unsigned  psizeO,u; 
float  slope, value; 

printf (" Select  mode  (0=manual , l=auto- enhance)  [0]:  ");  mode=getnum( ) 
if  (mode!=l)  mode=0; 

psizeO=psize [0] ; 
if  (mode==0)  { 
exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  lowest  value  (0,...,254):  ");  low=getnum( ) ; 
if  ( (low>=0)&&(low<=254) ) exit=l ; 

}; 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  highest  value  (%d, . . . , 255) : ”,low+l); 
high=getnum( ) ; 

if  ( (high>=(low+l) )&&(high<=255) ) exit=l ; 

); 

} ; 

if  (mode==l)  { 

low=255;  high=0; 
for  (u=0;  u<psize0;  ++u)  { 
if  (p[u]<low)  low=p[u]; 
if  (p[u]>high)  high=p [u] ; 

); 

printf ( "Lowest  value  is  %d,  Highest  value  is  %d . \n" , low , high) ; 

}; 

printf ("Enhancing  contrast. . . \n" ) ; 
slope=255 . 0/(high- low) ; 
for  (u=0;  u<psize0;  ++u)  { 
value=slope*(p [u] -low) ; 
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if  (value<0.0)  p[u]=0; 
else  { 

if  (value>255 . 0)  p[u]=255; 
else  p [u]=value+0 . 5 ; 

) ; 

); 

) 


fourier (p ,psize) 
char  *p; 

unsigned  Upsize; 

{ 

struct  complex  {float  real,imag;  }; 

struct  complex  *array , *getmem( ) , *getml ( ) , *temp ; 

char  line [MAXLINE] ; 

int  dir , exit , i , liml , lim2 , ret , scale ; 
unsigned  arraybytes , fftsize , beta , maxbeta , 

psizeO , psizel , psize2 , px , py , tempx, tempy ; 
long  tempbytes; 

float  beta4 , imag .highval , real , scale factor , val , root [ 256  ] ; 

printf("2-D  fast  fourier  transform  algorithm  (FFT)\n"); 
psizeO=psize [0] ; psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
tempbytes=8L*psize0 ; 
temp=getml( tempbytes) ; 

if  (temp==NULL)  printf ( "%c*****  insufficient  memory  *****\n",7); 
else  { 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  transform  direction  (f=foward, i=inverse)  [f]:  ") 
ret=getline(line) ; if  (ret==0)  line [ 0 ] = ' f ' ; 
switch  (line[0])  { 

case'F':  case'f':  dir=0;  exit=l ; break; 

case' I':  case'!':  dir=-l;  exit=l ; break; 

default:  printf ( "%c*****  invalid  entry  *****\n",7); 

} ; 

}; 

printf ( "Enter  power  of  2 fourier  plane  magnification  "); 
maxbeta=MAXFFT/psize2 ; 

printf ( " ( 1 , 2 , 4 , . . . , %d)  [1]:  ", maxbeta);  beta=getnum( ) ; 
ret=0;  for  (i=l;  i<=maxbeta;  i*=2)  if  (beta==i)  ret=l ; 
if  (ret==0)  beta=l ; 

f ftsize=psize2*beta;  beta4=( float ) bet a*beta*beta*be ta ; 
arraybytes=8*f ftsize ; array=getmem(arraybytes) ; 
if  (array==NULL)  printf ( "%c*****  malloc  error  *****\n" , 7 ) ; 
printf ("FFT  array  size  is  %d\n" , fftsize) ; 

printf ( "Calculating  square  root  table  ... \n" ) ; 
for  (i=0;  i<256;  ++i)  root [ i ]=sqrt( (double) i) ; 

printf ("Calculating  vertical  transforms. . .\n") ; 
printf("  0%%  completed\b\b\b\b\b\b\b\b\b\b\b\b\b\b " ) ; 
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liml=psizel/2+l ; 
lim2=psizel/2+(fftsize-psizel) ; 
tempx=0 ; 

while  (tempx<psize2)  { 

py=(psizel/2 -l)*psize2; 
for  (i=0;  i<fftsize;  ++i)  { 

if  ( ( i<liml ) | | ( i>lim2 ) ) array [ i ]. real=root [p [ tempx+py] ] ; 

else  array [ i] . real=0 . 0 ; 
array [ i ] . imag=0 . 0 ; 
py+=psize2;  if  (i==lim2)  py=0 ; 

}; 

fft87 (fftsize , dir , array) ; 
tempy=0 ; i=lim2+l ; 
while  ( tempy<psizeO)  { 

temp [ tempx+tempy ] . real=array [ i ] . real ; 
temp [ tempx+tempy ] . imag=array [ i ] . imag ; 
tempy+=psize2 ; 

++i ; if  (i==fftsize)  i=0 ; 

}; 

++tempx ; 

printf ( " %3 . Of\b\b\b" , ( 100 . 0*tempx) /psize2 ) ; 

); 

printf ( "\n" ) ; 

printf ( "Calculating  horizontal  transforms  ... \n" ) ; 
printf("  0%%  completed\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ) ; 
liml=psize2/2+l ; 
lim2=psize2/2+(fftsize-psize2) ; 
tempy=0 ; highval=0.0; 
while  ( tempy<psizeO)  { 
tempx=psize2/2 - 1 ; 
for  (i=0;  Kfftsize;  ++i)  { 
if  ( (icliml) | | (i>lim2) ) { 

array [i] . real=temp [ tempx+tempy ] .real; 
array [ i ] . imag=temp [ tempx+tempy ] . imag ; 

} 

else  { 

array [ i ] . real=0 . 0 ; array [ i ] . imag=0 . 0 ; 

); 

++tempx;  if  (i==lim2)  tempx=0 ; 

); 

fft87 (fftsize , dir , array) ; 
tempx=0 ; i=lim2+l ; 
while  ( tempx<psize2)  { 

real=array [ i ] . real ; imag=array [ i ] . imag ; 

val=(real*real+imag*imag)/beta4 ; temp [ tempx+tempy ] . real=val; 
if  (val>highval)  highval=val; 

++i;  if  (i==fftsize)  i=0 ; 

++tempx ; 

); 

tempy+=psize2 ; 

printf ( "%3 . 0f\b\b\b" , ( 100 . 0*tempy ) /psizeO) ; 

); 

printf ( "\n" ) ; 
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printf ( "Highest  value  is  % . 3e\n" ,highval) ; 

printf ( "Enter  integer  scale  magnification  ( l=autoscale)  [1] : ") 
scale=getnum( ) ; if  (scale<=0)  scale=l; 
printf ("Scaling. . .\n") ; 

printf ("  0%%  completed\b\b\b\b\b\b\b\b\b\b\b\b\b\b" ) ; 

scalefactor=255 . 0*scale/highval ; 
tempy=0 ; 

while  ( tempy<psizeO)  { 

for  (tempx=0;  tempx<psize2 ; -H-tempx)  { 

val=scalef actor*temp [ tempx+tempy ] . real+0 . 5 ; 
if  (val>255.0)  p [tempx+tempy ] =255 ; 
else  p [ tempx+tempy ]=val ; 

) ; 

tempy+=psize2 ; 

printf ( "%3 . Of\b\b\b" , (100 . 0*tempy)/psize0) ; 

}; 

printf ("\n") ; 

ret=rlsml ( temp , tempbytes) ; 
ret=rlsmem( array , arraybytes) ; 

}; 

) 

gaussian(p ,psize) 
char  *p; 

unsigned  *psize; 

{ 

char  c, line [MAXLINE] ; 
int  cutoff , exit , f ile , k, option, ret ; 
unsigned  index , loop, offset, psiz el, psize2,px,py,u; 
float  a , alpha, b , center , centerx , centery ,m, sumx, sumy , sumxy , sumx2 , 
value, x,y; 

long  sump , sumxp , sumyp ; 

printf ("Gaussian  Curve  Fit:\n"); 

exit=0 ; 

while  (exit=0)  { 

printf ( "Enter  maximum  pixel  noise  cutoff  value  (0,1,...)  [0]:  " 
cutof f=getnum( ) ; 

if  ( (cutof f>=0)&&(cutof f<=255) ) exit=l; 

else  printf  ("%c*****  invalid  value  *****\n"  , 7)  ; 

}; 

printf ( "Calculating  centroid. . .\n") ; 
sump=0 ; sumxp=0 ; sumyp=0 ; 
u=0 ; psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
for  (py=0;  py<psizel;  ++py)  { 
for  (px=0;  px<psize2;  ++px)  { 
c=p[u] ; ++u; 
if  (c>cutoff)  { 
sump+=c ; 
sumxp+=px'*c ; 
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sumyp+=py*c ; 


} ; 

) ; 

centerx=( (float) sumxp) /sump ; 
centery=( (float) sumyp) /sump ; 

printf ( "Centroid  is  ( % . 2f , % . 2f )\n" , centerx , centery ) ; 
exit=0 ; 

while  (exit==0)  { 

printf ( "Select  horizontal  or  vertical  fit  (h,v):  "); 
ret=getline(line) ; 
if  (ret ! =0)  { 

switch  (line[0])  { 
case 'H' : 

case'h':  loop=psize [ 2 ] ; offset=l ; option=3; 

center=centerx ; index=centery+0 . 5 ; 
exit=l ; break; 

case ' V' : 

case  V : loop=psize [ 1 ] ; of f set=psize [ 2 ] ; option=4; 
center=centery ; index=centerx+0 . 5 ; 
exit=l ; break; 

default:  printf ("%c*****  invalid  option  *****\n",7);  break 

) ; 

}; 

}; 


printf ("Fitting  to  a Gaussian. .. \n" ) ; 
f ile=creat ( "improc#2 . dat" , - 1 ) ; 

ret=write(f ile ,&option, 1) ; ret=write ( f ile ,&index, 1) ; 
c=loop-l;  ret=write(file,&c, 1) ; 
if  (option==3)  index*=psize [ 2 ] ; 

sumx=0 ; sumx2=0 ; sumy=0 ; sumxy=0 ; k=0 ; 
for  (u=0;  u<loop;  ++u)  { 

c=p[ (unsigned)u*offset+index] ; 
if  (c>cutoff)  { 

++k ; 

x=(u-center)*(u-center) ; 

y=log( (float)c) ; 

sumx+=x ; 

sumx2+=x*x ; 

sumy+=y ; 

sumxy+=x*y ; 

); 

}; 

m=(k*sumxy- sumx*sumy) /(k*sumx2-sumx*sumx) ; 
b=( sumx2*sumy - sumx*sumxy ) / (k*sumx2 - sumx*sumx) ; 
alpha=-m ; 
a=exp (b) ; 

printf("A  = %f  alpha  = %e  c = %f\n" , a , alpha , center ) ; 
printf ("Beam  spot  radius :\n"); 

printf ("\tl  over  e squared  = %f\n" , sqrt ( 2 . 0/alpha) ) ; 
printf ("\t5  percent  = %f\n" , sqrt(2 . 99573/alpha) ) ; 
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printf("\t4  percent  = %f\n" , sqrt (3 . 21888/alpha) ) ; 
printf("\t3  percent  = %f\n" , sqrt ( 3 . 50656/alpha) ) ; 
printf( "Press  RETURN  to  get  graph. \n");  ret=getline ( line) ; 

printf ( "Watch  graphics  monitor  ... \n" ) ; 
for  (u=0;  uCloop;  ++u)  { 

value=a*exp ( -alpha*(u-center)*(u-center) ) ; 
if  (value>255 . 0)  c=255;  else  c=value+0.5; 
if  (c=26)  c=25 ; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write (f ile , &c , 1 ) ; 

} ; 

for  (u=0;  uCloop;  ++u)  { 

c=p [ (unsigned) u*offset+index] ; 
if  (c— 26)  c=25 ; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write (f ile , &c , 1 ) ; 

); 

ret=close(file) ; 


getline ( s ) 
char  s [ ] ; 

{ 

int  c , i ; 

for  (i=0 ; i<MAXLINE- 1 &&  (c=getchar ( ) ) ! =E0F  &&  c!='\n';  ++i) 
s [ i ] =c ; 

s[i]='\0';  return(i); 


getnum( ) 

{ 

char  line [MAXLINE] ; 
int  i,n; 

i=getline (line) ; i=sscanf (line , "%d" , &n) ; 
if  (i==l)  return(n) ; else  return(O) ; 


graph_2d(p , psize) 
char  *p; 

unsigned  *psize; 

{ 

char  c , line [MAXLINE] ; 

int  exit , i , file , lim, option, ret ; 

unsigned  index , loop , of f , x ; 

printf ("2-D  line  intensity  graph\n"); 
exit=0 ; 

while  (exit==0)  { 

printf ( "Select  horizontal  or  vertical  scan  (h,v):  "); 
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ret=getline ( line) ; 
if  (ret !=0)  { 

switch  (line[0])  { 
case ' H ' : case ' h ' : 

lim=psize [ 1 ] - 1 ; loop=psize [ 2 ] ; off=l ; option=0; 
exit=l ; break; 
case ' V' : case ' v' : 

lim=psize [2] - 1 ; loop=psize [ 1 ] ; of f=psize [ 2 ] ; option=l ; 
exit=l ; break; 

default:  printf ( "%c*****  invalid  option  *****\n" , 7) ; break; 

); 

} ; 

) ; 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  index  value  (0 , 1 , . . . , %d) : " , lim) ; index=getnum( ) ; 
if  ( (index>=0)&&( index<=lim) ) exit=l ; 

else  printf ( "%c*****  illegal  value  *****\n" , 7 ) ; 

) ; 

printf ( "watch  graphics  monitor  ... \n" ) ; 
f ile=creat ( " improc#2 . dat" , - 1 ) ; 

ret=write ( f ile , &option, 1 ) ; ret=write ( f ile , &index , 1 ) ; 
c=loop-l;  ret=write(file,&c, 1) ; 
if  (option==0)  index*=psize [ 2 ] ; 
for  (x=0;  x<loop;  ++x)  { 

c=p [ (unsigned) x*off+index] ; 
if  (c— 26)  c=25  ; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write (f ile , &c , 1 ) ; 

}; 

ret=close(file) ; 

} 

histogram(p ,psize) 
char  *p ; 

unsigned  upsize; 

{ 

int  ave , f ile , i , j , highval , lowval , option , peakval , ret ; 
unsigned  u , peak ,h256 [ 256 ] ; 
long  sum; 

printf ( "Computing  histogram. . . \n") ; 
for  (i=0;  i<256;  ++i)  h256[i]=0; 
for  (u=0;  u<psize[0];  ++u)  ++h256 [p [u] ] ; 
sum=0 ; 

for  (i=0;  i<2 5 6 ; ++i)  sum+=( long) i*h256  [ i] ; 
peak=0 ; 

for  (i=0;  i<2 5 6 ; ++i)  if  (h256 [ i]>peak)  { peak=h256 [ i ] ; peakval=i ; ); 
highval=255 ; 

while  (h256 [highval ] ==0)  --highval; 
lowval=0 ; 

while  (h256 [ lowval ] ==0)  ++lowval; 
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ave-(int) (( float ) sum/psize [ 0 ] ) ; 

printf ("Creating  histogram  graph  (watch  graphics  monitor )... \n" ) ; 
f ile=creat ( " improc#2 . dat" , - 1 ) ; 
option=2 ; ret=write(file ,&option, 1) ; 

ret=write(file ,&lowval , 1) ; ret=write(file ,&highval , 1) ; 
ret=write(f ile , &ave , 1 ) ; ret=write ( f ile , &peakval , 1 ) ; 
for  (i=0;  i<256;  ++i)  { 

j = ( ( long) 255*h256 [ i ] ) /peak ; 
if  (j==0  &&  h256 [ i ] >0)  j=l; 
if  ( j ==26 ) j=25 ; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write(file,&j , 1) ; 

) ; 

ret=close(file) ; 


init_pic (p , psize) 
char  *p ; 

unsigned  *psize; 

{ 

int  i, file, ret; 

f ile=open( " improc#l . pic" , 0_RD0NLY  | 0x8000); 
if  (f ile==- 1 ) { 

/*  interupted  picture  file  does  not  exist,  begin  fresh  */ 
printf ("\nIMAGE  PROCESSING  PR0GRAM\n" ) ; 
printf("by  Matt  Weppner  (NBS  1986)\n"); 
for  (i=0;  KPSIZESIZE;  ++i)  psize[i]=0; 

) 

else  { 

/*  interupted  picture  file  is  retrieved  */ 
psize [0]=0;  ret=FBYTES; 

while  (psize [0]<MAXPSIZE  &&  ret==FBYTES)  { 

ret=read(f ile ,p , FBYTES) ; p+=ret;  psize [ 0 ] +=ret ; 

); 

calc_res (psize) ; 

); 

ret=close(f ile) ; ret=unlink( " improc#l .pic") ; 

} 

input_matrix(psize ,ma,mv) 
unsigned  *psize; 
int  *ma,*mv; 

{ 

int  a, i , j ,n,v,x,y; 

printf ( "Convolution  with  input  matrix\n"); 

printf ( "Enter  matrix  height  (1,3,...):  ");  y=getnum(); 

printf ( "Enter  matrix  width  (1,3,...):  ");  x=getnum(); 

printf ( "\nRow\tColumn\tValue\n" ) ; 

n=0 ; 

for  (j=l ; j <=y ; ++j  ) { 
for  (i=l;  i<=x;  ++i)  { 
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printf ( " %d\t%d\t" , j , i) ; v=getnum( ) ; 
if  (v!=0)  { 

++n; 

ma[n] =psize [ 2 ] *( j - (y/2) - 1 ) + ( i - (x/2) - 1 ) ; mv[n]=v; 


printf("\n%d  non-zero  elements . \n" , n) ; ma[0]=n; 

} 


mag(pp,psize) 
char  **pp ; 
unsigned  Upsize; 

{ 

char  c , *malloc ( ) , *p , *temp ; 
int  ret; 

unsigned  psizel ,psize2,px, pxstart , pxf inish , py , py start , py finish , 
tempx , tempy ; 

printf ( "Magnifying  by  2X. . . \n") ; 
p=*pp ; temp=malloc (MAXPSIZE) ; 

if  (temp==NULL)  printf (" %c*****  malloc  error  *****\n",7); 
else  { 

psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
pxstart=psize2/4 ; pxf inish=(3*psize2) /4 ; 

pystart=(psize 1/4) Upsize 2 ; pyfinish=( (3*psizel) /4)*psize2 ; 
tempy=0 ; 

for  (py=pystart;  py<pyf inish;  py+=psize2)  { 
tempx=0 ; 

for  (px=pxstart;  px<pxf inish;  ++px)  { 
c=p[px+py] ; 

temp [ tempx+tempy ] =c ; temp [ tempx+l+tempy ] =c ; 

temp [ tempx+tempy+psize2 ] =c ; temp [ tempx+l+tempy+psize2 ] =c ; 

tempx+=2 ; 

); 

tempy+=(2*psize2) ; 

); 

*pp=temp;  ret=free(p); 

} ; 

} 


multi_mode(p ,psize) 
char  *p ; 

unsigned  Upsize; 

{ 

char  c , 1 ine [ MAXLINE ] , pmax , pmin ; 

int  cutoff , exit, file, iO, imax , k , n, option , ret , radius , step ; 
unsigned  index , loop, offset, psizel ,psize2 ,px,py,u; 

float  a ,b , center , centerx , centery , g , e , m , r , sumx , sumy , sumxy , sumx2 , value , x , y 
float  olda , olde , oldg , rarray [ 256 ] ,xarray[256] ,yarray[256] ; 
long  sump , sumxp , sumy p ; 

printf ( "Multi -mode  Curve  Fit:\n"); 
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exit-=0 ; 

while  (exit==0)  { 

printf( "Enter  maximum  pixel  noise  cutoff  value  (0,1,...)  [0]:  "); 
cutof f=getnum( ) ; 

if  ( (cutoff>=0)&&(cutoff<=255) ) exit=l ; 

else  printf ("%c*****  invalid  value  *****\n",7); 

} ; 

printf ( "Calculating  centroid. . . \n") ; 
sump=0 ; sumxp=0 ; sumyp=0 ; 
u=0;  psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
for  (py=0;  py<psizel;  ++py)  { 
for  (px=0;  px<psize2;  ++px)  { 
c=p[u] ; ++u; 
if  (c>cutoff)  { 
sump+=c ; 
sumxp+=px*c ; 
sumyp+=py*c ; 

) ; 


centerx=( ( float )sumxp) /sump ; 
centery=( (float) sumyp) /sump ; 

printf ( "Centroid  is  (% . 2f , % . 2f)\n" , centerx, centery) ; 
exit=0 ; 

while  (exit==0)  { 

printf ("Select  horizontal  or  vertical  fit  (h,v):  "); 
ret=getline ( line) ; 
if  (ret ! =0)  { 

switch  (line [0] ) { 
case ' H ' : 

case'h':  loop=psize [ 2 ] ; offset=l ; option=5; 

center=centerx ; index=centery+0 . 5 ; 
exit=l ; break; 

case ' V' : 

case'v' : loop=psize [ 1 ] ; offset=psize [ 2 ] ; option=6; 
center=centery ; index=centerx+0 . 5 ; 
exit=l;  break; 

default:  printf ("%c*****  invalid  option  *****\n",7);  break; 

) ; 

); 

}; 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  radius  of  exclusion  (0,1,...)  [0]:  ");  radius=getnum( ) 
if  ( (radius>=0)&&(radius<=loop/3) ) exit=l ; 
else  printf ( "%c*****  invalid  value  *****\n" , 7 ) ; 

); 

printf ( "Fitting  to  a g-prof ile . . . \n" ) ; 
f ile=creat ( "improc#2 . dat" , - 1 ) ; 
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ret=write(file ,&option, 1) ; ret=write ( f ile ,&index, 1) ; 
c=loop-l;  ret=write(file ,&c , 1) ; 
if  (option==5)  index*=psize [ 2 ] ; 


imax=0 ; 

for  (u=0;  u<loop;  ++u)  { 

c=p [ (unsigned)u*of fset+index] ; 
if  (c>imax)  imax=c ; 

} ; 

printf("imax  = %d\n" , imax) ; 

pmax=imax*EIAMMMAX ; pmin=imax*EIAMMMIN ; 
k=0 ; sumx=0 ; sumx2=0 ; 
for  (u=0;  u<loop ; ++u)  { 

c=p [ (unsigned)u*offset+index] ; r=fabs (u- center) ; 
if  ( (c>=pmin)&&(c<=pmax)&&(r>radius) ) { 

rarray[k]=r;  x=log(r) ; xarray[k]=x;  yarray[k]=c; 
sumx+=x  ; sumx2+=x*x;  ++k; 

} ; 

) ; 

printf("%d  data  points  used.\n",k); 

printf ( "Enter  starting  value  for  Io  [ % d ] : " , imax) ; iO=getnum(); 
if  ( (i0==0) | | (i0>=255) ) iO=imax; 
n=0 ; iO+=l ; step=-l; 
e=1.0el0;  a=0 ; g=0 ; 

while  ( ( (e<olde)  | | (n<=2)  )&&(i0>=0)&&(i0<=255) ) { 
olde=e;  olda=a;  oldg=g; 

++n;  iO+=step; 
sumy=0 ; sumxy=0 ; 
for  (u=0;  u<k;  ++u)  { 

x=xarray[u] ; y=log(l . 0-yarray [u] /iO) ; 
sumy+=y;  sumxy+=x*y; 

}; 

m=(k*sumxy- sumx*sumy )/(k*sumx2 - sumx*sumx) ; g=m; 
b=(sumx2*sumy-sumx*sumxy)/(k*sumx2-sumx*sumx) ; a=exp(-b/m) ; 
e=0  ; 

for  (u=0;  u<k;  ++u)  { 

value=yarray [u] -iO*(l . 0 - exp (g*log(rarray [u] /a) ) ) ; 
e+=value*value ; 

) ; 

printf ("n=%d  Io=%d  error=% . 4e\n" , n , iO , e) ; 
if  (e>olde)  step=-step; 

} ; 

iO+=step;  a=olda;  g=oldg; 

printf ( "\nFinal  Parameters:  Io  = %d  a = %f  (pixels)  g = %f\n\n" , iO , a , g) 
printf ( "Press  RETURN  to  get  graph. \n");  ret=getline ( line) ; 

printf ( "Watch  graphics  monitor  ... \n" ) ; 
for  (u=0;  u<loop;  ++u)  { 

value=iO* ( 1 . 0 -exp ( g*log( fabs (u- center ) /a) ) ) ; 
if  (value>255 . 0)  c=255; 

else  { if  (value<0.0)  c=0 ; else  c=value+0 . 5 ; }; 

if  (c==26)  c=25 ; 
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/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write(file , &c , 1) ; 

); 

for  (u=0;  u<loop ; ++u)  { 

c=p[ (unsigned)u*of fset+ index] ; 
if  ( c==26 ) c=25; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC  which 
halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
ret=write(file ,&c , 1) ; 

) ; 

ret=close(file) ; 


negative (p , psize) 
char  *p; 

unsigned  *psize; 

{ 

unsigned  psizeO,u; 

printf( "Creating  negative  picture  (reversing  contrast) ... \n" ) ; 
psizeO=psize [0] ; 

for  (u=0;  u<psize0;  ++u)  p [u]=255-p [u] ; 


plot_3d(p , psize) 
char  *p ; 

unsigned  *psize; 

{ 

char  c , lowest , highest , line [MAXLINE] , *temp , *malloc ( ) ; 
int  exit , file , i , ret , zoom; 
unsigned  center , tempsize , t ,u,v, 

x , y ,xsize , ysize , xcenter , ycenter , xmin , ymin , xmax , ymax , yj  ump ; 

printf("3-D  intensity  graph\n"); 

zoom=0 ; 

if  (psize [ 3]>=4)  { 
exit=0 ; 

while  (exit==0)  { 

printf( "Select  full  size  or  4X  enlargement  (f,e):  "); 
re  t=ge  1 1 ine ( 1 ine ) ; 
if  (ret ! =0)  { 

switch  (line[0])  { 
case ' F'  : 

case'f':  exit=l ; break; 
case ' E ' : 

case'e' : zoom=l ; exit=l ; break; 

default:  printf ("%c*****  invalid  entry  *****",7);  break: 

); 

}; 

) ; 

); 
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if  (zoom==l)  { 

xmin=psize [ 2 ] /8 - 1 ; xmax=psize [ 2 ] - (xmin+2) ; center=psize [ 2 ] /2 - 1 ; 
printf( "Enter  center  horizontal  position  "); 
printf ( " ( %d , . . . , %d)  [%d]:  " ,xmin,xmax, center) ; 
xcenter=getnum( ) ; 

if  ( (xcenter<xmin) | | (xcenter>xmax) ) xcenter=center ; 
ymin=psize [ 1 ] /8 - 1 ; ymax=psize [ 1 ] - (ymin+2 ) ; center=psize [ 1 ] /2 - 1 ; 
printf ( "Enter  center  vertical  position  "); 
printf (" (%d %d)  [%d]:  ", ymin , ymax , center ) ; 

ycenter=getnum( ) ; 

if  ( (ycenter<ymin) | | (ycenter>ymax) ) ycenter=center ; 

} ; 

printf ( "Creating  3-D  plot  (watch  graphics  monitor) ... \n" ) ; 
switch  ( ( int)psize [ 3 ] ) { 


case 

1 : 

ysize=15 

xsize=16 ; 

break; 

case 

2: 

ysize=30 

xsize=32 ; 

break ; 

case 

3: 

ysize=60 

xsize=64 ; 

break; 

case 

4: 

if  (zoom= 

==0)  {ysize= 

=60;  xsize=128;  } 

else  { ysize=30;  xsize=32;  }; 
break; 


case  5:  ysize=60; 

if  (zoom==0)  { xsize=256;  } else  { xsize=64;  }; 
break; 

} ; 

yjump=(psize [ 1 ] /ysize) Upsize [ 2 ] ; 

tempsize=ysize*xsize ; temp=malloc( tempsize) ; 

if  (temp==NULL)  printf ( "*****  malloc  error  *****\n"); 

f ile=creat ( " improc#3 . dat" , -1) ; 

c=ysize-l;  ret=write ( f ile , &c , 1 ) ; 

c=xsize-l;  ret=write (f ile , &c , 1)  ; 

ret=write (f ile , &zoom, 1 ) ; 

lowest=255;  highest=0; 

if  (zoom==0)  { 

t=0 ; u=psize [ 0 ] -psize [ 2 ] ; 
for  (y=0;  y<ysize;  ++y)  { 
for  (x=0;  x<xsize;  ++x)  { 

c=p[u+x];  if  (c==26)  c=25;  temp[t]=c;  ++t; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC 
which  halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
if  (cClowest)  lowest=c;  if  (c>highest)  highest=c; 

) ; 

u-=yj ump; 

); 

ret=write(f ile ,&lowest , 1) ; ret=write(file,&highest, 1) ; 
ret=write ( f ile , temp , tempsize) ; 

} 

else  ( 

t=0 ; u=psize [ 2 ] *(ycenter+psize [ 1 ]/8) ; xcenter-=psize [ 2 ] /8 - 1 ; 
for  (i=0;  i<ysize;  ++i)  { 
for  (v=0;  v<xsize;  ++v)  { 

c=p [u+v+xcenter ] ; if  (c==26)  c=25;  temp[t+v]=c; 

/*  Note:  this  step  is  taken  because  character  26  is  EOF  in  BASIC 
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which  halts  the  reading  of  a file  (ie.  26  is  illegal)  */ 
if  (cClowest)  lowes t=c ; if  (c>highest)  highest=c; 

); 

t+=xsize;  u-=psize[2]; 

) ; 

ret=wr ite ( f ile , &lowes t , 1) ; ret=write(file , &highest , 1) ; 
ret=write (f ile , temp , tempsize) ; 

) ; 

ret=close (f ile) ; ret=f ree ( temp) ; 


pratt (pp , psize , ma ,mv) 
char  **pp ; 
unsigned  *psize; 
int  *ma,*mv; 

{ 

int  i, epsilon; 

printf ( "Pratt ' s noise  cleaning  algorithim\nEnter  epilson:  "); 
epsilon=getnum( ) ; 
ma [ 0 ] =8 ; 
i=psize [ 2 ] ; 

ma [ 1 ] =- ( i+1 ) ; ma[2]=-i;  ma[3]=- (i-1) ; ma[4]=-l; 
ma[5]=l;  ma[6]=i-l;  ma[7]=i;  ma[8]=i+l; 
for  (i=l;  i<=8 ; ++i)  mv[i]=l; 
convolve (pp , psize ,ma,mv, 2 , epsilon) ; 


read_pic (p , psize) 
char  *p ; 

unsigned  *psize; 

{ 

char  line [MAXLINE] ; 
int  file, ret; 

printf ("Enter  input  filename:  ");  ret=getline(line) ; 
printf ("Reading  file  %s\n",line); 
f ile=open( line , 0_RD0NLY  | 0x8000); 

if  (file=-l)  printf ("%c*****  cannot  open  file  *****\n" , 7) ; 
else  { 

psize [ 0 ] =0 ; ret=FBYTES ; 

while  (psize [0]<MAXPSIZE  &&  ret==FBYTES)  { 

ret=read(f ile , p , FBYTES) ; p+=ret;  psize [ 0 ] +=ret ; 

}; 

printf ("%u  bytes  read . \n" , psize [ 0 ]) ; 
calc_res (psize) ; 

printf ( "Resolution  is  %ux%u\n" , psize [ 2 ], psize [ 1 ]) ; 

}; 

ret=close(file) ; 

} 

resolution(pp , psize) 
char  **pp ; 
unsigned  *psize; 
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char  *malloc ( ) , *p , *temp ; 
int  exit,i,ret; 

unsigned  area , ratio , sum , tempsize [ PSIZESIZE  ] , u , x , xof f , y , yof f , 
psizeO,psizel,psize2,tsizeO,tsizel,tsize2; 

printf( "Resolution  is  %ux%u\n" , psize [ 2 ] , psize [ 1 ] ) ; 
exit=0 ; 

while  (exit==0)  { 

printf( "Enter  new  resolution  "); 

printf ( " ( 1=16x15 , 2=32x30 , 3=64x60 , 4=128x120 , 5=256x240) : ") ; 
ret=getnum( ) ; 
switch  (ret)  { 


case  1 

tempsize [0]=  240 

exit=l 

break; 

case  2 

tempsize[0]=  960 

exit=l 

break ; 

case  3 

tempsize[0]=  3840 

exit=l 

break; 

case  4 

tempsize [0]=15360 

exit=l 

break; 

case  5 

tempsize [ 0 ] =6 1 440 

exit=l 

break; 

default:  printf  ("  ifec*****  invalid  option  **■***",  7);  break; 

}; 

) ; 

calc_res (tempsize) ; 

if  ( tempsize [ 3 ]! =psize [ 3 ] ) { 

p=*pp ; temp=malloc (MAXPSIZE) ; 

if  (temp==NULL)  printf ( "ifec*****  malloc  error  *****\n" , 7 ) ; 
psizeO=psize [ 0 ] ; psizel=psize [ 1 ] ; psize2=psize [ 2 ] ; 
tsizeO=tempsize [ 0 ] ; tsizel=tempsize [ 1 ] ; tsize2=tempsize [ 2 ] ; 
printf ( "Changing  resolution  from  %ux%u  ", psize [ 2 ], psize [ 1 ]) ; 
printf ( "to  %ux%u. . . \n" , tempsize [ 2 ] , tempsize [ 1 ] ) ; 

if  ( tsizeO<psizeO)  { 
ratio=psizel/tsizel ; 
area=ratio*ratio ; 
y=0 ; x=0 ; 

for  (u=0;  u<tsize0;  ++u)  { 
sum=0 ; 

for  (yoff=0;  yoff<(ratio*psize2) ; yoff+=psize2)  { 
for  (xoff=0;  xoff<ratio;  ++xoff)  { 
sum+=p [y+yof f+x+xoff ] ; 

); 

} ; 

temp [u] =sum/area ; 
x+=ratio ; 

if  (x==psize2)  { x=0 ; y+=ratio*psize2 ; ); 

} ; 

} ; 

if  ( tsizeO>psizeO)  { 
ratio=tsizel/psizel ; 
y=0;  x=0 ; 

for  (u=0;  u<psize0;  ++u)  { 

for  (yoff=0;  yoff<(ratio*tsize2) ; yoff+=tsize2)  { 
for  (xoff=0;  xoff<ratio;  ++xoff)  { 
temp [y+yof f+x+xoff ] =p [u] ; 

) ; 
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} ; 

x+=ratio ; 

if  (x==tsize2)  (x=0;  y+=ratio*tsize2 ; }; 

} ; 

) ; 

*pp=temp ; ret=free(p); 

for  (i=0;  i<PSIZESIZE;  ++i)  psize [ i ] =tempsize [ i ] ; 

) ; 

printf ( "Resolution  is  %ux%u\n" , psize [ 2 ], psize [ 1 ]) ; 


savepic (p , psize) 
char  *p; 

unsigned  *psize; 

{ 

int  file, ret; 
unsigned  sum; 

file=creat("improc#l .pic" , -1) ; 
sum=0 ; ret=FBYTES; 

while  (sum<psize [ 0 ] &&  ret==FBYTES)  { 

ret=write(file ,p , FBYTES) ; p+=ret;  sum+=ret; 

) ; 

ret=close(file) ; 

) 

threshold(p , psize) 
char  *p ; 

unsigned  *psize; 

{ 

int  exit , low , high , new; 
unsigned  psizeO.u; 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  lower  bound  [0, . . . ,255] : ");  low=getnum( ) ; 
if  ( ( low>=0)&&( low<=255) ) exit=l ; 

}; 

exit=0 ; 

while  (exit==0)  { 

printf ( "Enter  upper  bound  [ %d, . . . , 255 ] : " , low) ; high=getnum( ) 
if  ( (high>=low)&6c(high<=255)  ) exit=l ; 

) ; 

exit=0 ; 

while  (exit==0)  { 

printf ("Enter  new  value  [0,...,255]:  ");  new=getnum( ) ; 
if  ( (new>=0)&&(new<=255) ) exit=l ; 

); 

printf ( "Setting  threshold  values  ... \n" ) ; 

psizeO=psize [0 ] ; 

for  (u=0;  u<=psize0;  ++u) 

if  ( (p [u]>=low)&&(p [u]<=high) ) p[u]=new; 
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write_pic (p , psize) 
char  *p ; 

unsigned  Upsize; 

{ 

char  Line [MAXLINE] ; 
int  file, ret; 
unsigned  sum; 

printf( "Enter  output  filename:  ");  ret=getline ( line ) ; 
printf( "Writing  file  %s\n",line); 
f ile=creat ( line , - 1 ) ; 

if  (file==-l)  printf ( " %c*****  cannot  open  f ile*****\n" , 7 ) 
else  { 

sum=0 ; ret=FBYTES; 

while  (sum<psize [0]  &&  ret==FBYTES)  { 

re t=write ( f ile , p , FBYTES) ; p+=ret;  sum+=ret; 

} ; 

} ; 

printf ("%u  bytes  written . \n" , sum) ; 
ret=close (f ile) ; 
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IMPR0C#2 . BAS 


10  SCREEN  2 : CLS 
20  KEY  OFF 

30  OPEN  " improc#2 . dat " FOR  INPUT  AS  1 
40  XSCALE=2 : XB=76 
50  YSCALE= . 5 : YB=20 
60  0PT=ASC ( INPUT$ ( 1 , #1 ) ) 

70  PRINT  "2-D  PLOTTING  PROGRAM" 

80  PRINT 

90  PRINT  "When  plot  is  completed  press  Shift-Prt  to  print  hardcopy," 
100  PRINT  "or  any  other  key  to  exit." 

110  PRINT 
120  MODE=0 

130  IF  ( OPT>=3 ) AND ( OPT<=6 ) GOTO  160 

140  INPUT  "Select  Graph  Mode  (0=line , l=bar) [0] : " .MODE 
150  PRINT 

160  IF  MODEOl  THEN  MODE=0 

170  IF  (MODE=0 ) OR ( 0PTO2 ) THEN  GOTO  230 

180  INPUT  "Enter  Bar  Width  (odd  # of  pixels  only)  [ 1 ] : " , BARWIDTH 
190  PRINT 

200  BARWIDTH=CINT( BARWIDTH) 

210  IF  BARWIDTHCl  THEN  BARWIDTH=1 
220  BARWIDTH=INT(BARWIDTH/2) 

230  INPUT  "Message: " ,MESSAGE$ 

240  CLS 

250  IF  0PTO0  THEN  GOTO  270 

260  PRINT "Horizontal  Intensity  Graph:  " ; :XLIM=255 : GOTO  450 
270  IF  OPTOl  THEN  GOTO  290 

280  PRINT "Vertical  Intensity  Graph:  " ; :XLIM=239 : GOTO  450 
290  IF  0PTO2  THEN  GOTO  360 

300  L0W=ASC ( INPUT $ ( 1 , #1 ) ) : HIGH=ASC ( INPUT? ( 1 , #1 ) ) 

310  AVE=AS C( INPUT? ( 1,#1) ) :PEAK=ASC( INPUT? (1,#1)) 

320  XLIM=255 : LOOP=255 

330  PRINT  "Intensity  Histogram:  Low=" ; LOW; " High=";HIGH; 

340  PRINT  " Average=" ;AVE; " Peak=";PEAK 

350  GOTO  480 

360  IF  0PTO3  THEN  380 

370  PRINT"Horizontal  Gaussian  Curve  Fit:  " ; :XLIM=255 : GOTO  450 
380  IF  0PTO4  THEN  GOTO  400 

390  PRINT "Vertical  Gaussian  Curve  Fit:  " ; :XLIM=239 : GOTO  450 
400  IF  0PTO5  THEN  GOTO  420 

410  PRINT "Horizontal  G-Profile  Fit:  " ; :XLIM=255 : GOTO  450 
420  IF  0PTO6  THEN  GOTO  440 

430  PRINT "Vertical  G-Profile  Fit:  " ; :XLIM=239 : GOTO  450 
440  PRINT  "Illegal  input  file": STOP 
450  INDEX=ASC (INPUT? ( 1 ,#1)) 

460  LOOP=ASC( INPUT? (1 ,#1) ) 

470  PRINT  "Line  Index  ="; INDEX 
480  PRINT 

490  IF  OPT=2  THEN  GOTO  550 
500  FOR  Y=256  TO  0 STEP  -64 
510  PRINT  TAB (6) ;Y 

520  PRINT : PRINT : PRINT 
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530  NEXT  Y 
540  GOTO  600 

550  PRINT  TAB(5)"1 .00" : PRINT: PRINT: PRINT 
560  PRINT  TAB (5)" 0.7 5": PRINT: PRINT: PRINT 
570  PRINT  TAB (5)" 0.50": PRINT: PRINT: PRINT 
580  PRINT  TAB ( 5) "0. 25" : PRINT: PRINT: PRINT 
590  PRINT  TAB (5) "0 . 00" 

600  LINE  (XB , YB) - (XB+2*XLIM+2 , YB) 

610  LINE  (XB+2*XLIM+2,YB)-(XB+2*XLIM+2,YB+128) 

620  LINE  (XB , 128+YB) - (XB+2*XLIM+2 , 128+YB) 

630  LINE  (XB , YB) - (XB , YB+128) 

640  SSIZE=16 :TSIZE=1 

650  FOR  X=0  TO  (2*XLIM)+2  STEP  SSIZE 

660  LINE(XB+X, YB) - (XB+X, YB-TSIZE) 

670  LINE (XB+X, YB+128) - (XB+X , YB+128+TSIZE) 

680  NEXT  X 

690  SSIZE=8 : TSIZE=3 

700  FOR  Y=0  TO  128  STEP  SSIZE 

710  LINE(XB, YB+128 -Y) - (XB-TSIZE , YB+128 -Y) 

720  LINE (XB+2*XLIM+2 , YB+1 2 8 - Y) - (XB+2*XLIM+2+TS IZE , YB+1 2 8 - Y) 

730  NEXT  Y 

740  PRINT  CHR$ (11) 

750  IF  (OPT=0)OR(OPT=1)OR( (OPT>=3)AND(OPT<=6) ) THEN  Y$=" INTENSITY" : Zl=5 : Z2=4 
760  IF  (OPT=2)  THEN  Y$="REL  FREQUENCY" : Z 1=3 : Z2=2 
770  FOR  Y=1  TO  Zl : PRINT : NEXT  Y 

780  FOR  Y=1  TO  LEN(Y$ ): PRINT  MID$ (Y$ , Y , 1 ) : NEXT  Y 
790  FOR  Y=1  TO  Z2 : PRINT : NEXT  Y 


800 

IF 

(LOOP=255) 

THEN  PRINT  TAB (10) 

II 

0 

64 

128 

192 

256" 

810 

IF 

(LOOP=127) 

THEN  PRINT  TAB (10) 

!t 

0 

32 

64 

96 

128" 

820 

IF 

(LOOP=63) 

THEN  PRINT  TAB (10) 

It 

0 

16 

32 

48 

64" 

830 

IF 

(LOOP=31) 

THEN  PRINT  TAB (10) 

II 

0 

8 

16 

24 

32" 

840 

IF 

(LOOP=15) 

THEN  PRINT  TAB (10) 

II 

0 

4 

8 

12 

16" 

850 

IF 

(LOOP=239) 

THEN  PRINT  TAB (10) 

II 

0 

40  80 

120 

160 

200 

240" 

; 

860 

IF 

(LOOP=119) 

THEN  PRINT  TAB (10) 

It 

0 

20  40 

60 

80 

100 

120" 

} 

870 

IF 

(LOOP=59) 

THEN  PRINT  TAB (10) 

II 

0 

10  20 

30 

40 

50 

60" 

; 

880 

IF 

(LOOP=29) 

THEN  PRINT  TAB (10) 

II 

0 

5 10 

15 

20 

25 

30" 

i 

890 

IF 

(L00P=14) 

THEN  PRINT  TAB (10) 

II 

0 

5 

10 

15" 

> 

900  PRINT 

910  IF  (OPT=0)OR(OPT=3)OR(OPT=5)  THEN  PRINT  TAB(35);"P  0 S I T I 0 N" 
920  IF  (OPT=l ) OR(OPT=4)OR(OPT=6)  THEN  PRINT  TAB(33);"P  0 S I T I 0 N" 
930  IF  (OPT=2)  THEN  PRINT  TAB(34);"I  N T E N S I T Y" 

940  PRINT: PRINT  MESSAGE? ; 

950  FIRST=1 

960  IF  (OPT=Q)OR(OPT=3)OR(OPT=5)  THEN  XSCALE=512/(LOOP+l ) 
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970  IF  (OPT=l )OR(OPT=4)OR(OPT=6)  THEN  XSCALE=480/ (LOOP+1 ) 

980  FOR  X=0  TO  LOOP 

990  Y=ASC ( INPUT$ ( 1 , #1 ) ) *YSCALE 

1000  XSCREEN=XSCALE*X+XB : YSCREEN=128 -Y+YB- . 1 

1010  IF  MODE=l  THEN  GOTO  1040 

1020  IF  FIRST=0  THEN  LINE  - (XSCREEN , YSCREEN) : GOTO  1080 
1030  PSET(XSCREEN, YSCREEN) :FIRST=0: GOTO  1080 

1040  IF  (0PTO2)  THEN  LINE  (XSCREEN , YB+128 )-  (XSCREEN , YSCREEN)  : GOTO  1080 
1050  FOR  I=-BARWIDTH  TO  BARWIDTH 

1060  LINE(XSCREEN+I , YB+128) - (XSCREEN+I , YSCREEN) 

1070  NEXT  I 

1080  NEXT  X 

1090  IF  NOT ( ( OPT>=3 ) AND ( OPT<=6 ) ) THEN  GOTO  1220 

1100  1=1 

1110  FIRST=1 

1120  FOR  X=0  TO  LOOP 

1130  Y=ASC(INPUT$(1 ,#1))*YSCALE 

1140  XSCREEN=XSCALE*X+XB:YSCREEN=128-Y+YB- .1 

1150  IF  FIRST=1  THEN  FIRST=0 : PSET (XSCREEN , YSCREEN) : GOTO  1210 

1160  IF  (LOOPMOO)  THEN  LINE- (XSCREEN  .YSCREEN)  : GOTO  1210 

1170  LINE (XSCREEN - 2*1 , YSCREEN+I ) - (XSCREEN+2*I , YSCREEN+I ) 

1180  LINE(XSCREEN-2*I .YSCREEN-I) - (XSCREEN+2*I , YSCREEN- I ) 

1190  LINE (XSCREEN+2*I , YSCREEN-I ) - (XSCREEN+2*I , YSCREEN+I ) 

1200  LINE(XSCREEN-2*I , YSCREEN-I) - (XSCREEN- 2*1 , YSCREEN+I) 

1210  NEXT  X 

1220  CLOSE  1 

1221  BEEP 

1230  A$=INKEY$ : IF  A$=""  THEN  1230 
1240  KILL  " improc#2 . dat" 

1250  SYSTEM 
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IMPR0C#3 . BAS 


10  SCREEN  2 : CLS 
20  KEY  OFF 

30  PRINT  "3-D  PLOTTING  PROGRAM" 

40  PRINT 

50  PRINT  "When  plot  is  completed  press  Shift-PrtSc  to  print  a hardcopy," 

60  PRINT  "or  any  other  key  to  exit." 

70  PRINT 

80  INPUT  "Select  Autoscale  (l=on,0=off)  [1]:  ",AUTOSCALE$ 

90  PRINT 

100  IF  AUTOSCALE$="0"  THEN  AUTOSCALE=0  ELSE  AUT0SCALE=1 

110  INPUT  "Select  Graph  Type  (0=lines , l=cross -hatch)  [0]:  ",TYPE$ 

120  PRINT 

130  IF  TYPE$=" 1 " THEN  TYPE=1  ELSE  TYPE=0 
140  INPUT  "Message:  " , MESSAGES 
150  CLS 

160  PRINT  "THREE  DIMENSIONAL  INTENSITY  PLOT:  " ; MESSAGE$ 

170  DIM  BIGGEST ( 639 ) , SMALLEST (6 39) , TEMPBIGGEST(639) , TEMPSMALLEST (639) 

180  OPEN  " improc#3 . dat"  FOR  INPUT  AS  1 

190  YSIZE=ASC (INPUT$ ( 1 , #1 ) ) :XSIZE=ASC(INPUT$(1 ,#1) ) 

200  DIM  OLDLINE (256,2) 

210  ZOOM=ASC (INPUT $ ( 1 , #1 ) ) 

220  LOWEST=ASC ( INPUT$ ( 1 , #1 ) ) : HIGHEST=ASC ( INPUT$ ( 1 , #1 ) ) 

230  IF  AUT0SCALE=0  THEN  LOWEST=0 : HIGHEST=255 

240  MSCALE=255/ (HIGHEST-LOWEST) : BSCALE=(255*LOWEST)/(LOWEST-HIGHEST) 

250  HEIGHT= . 27 : FIRSTX=0 : LASTX=XSIZE 
260  FOR  ELEMENTS  TO  639 

270  BIGGEST ( ELEMENT )=0 : SMALLEST (ELEMENT) =19 9 

280  TEMPBIGGEST ( ELEMENT ) =0 : TEMPSMALLEST ( ELEMENT) =1 9 9 

290  NEXT  ELEMENT 

300  XSTEP=2*(256/(XSIZE+1)) 

310  YSTEP=2*(60/(YSIZE+1)) 

320  FOR  Y=0  TO  YSIZE 
330  XSCREENADJUST=6+YSTEP*Y 

340  YSCREENADJUST=199-YSTEP*Y 

350  FOR  X=FIRSTX  TO  LASTX 

360  Z= (HEIGHT* (MSCALE*ASC ( INPUT$ ( 1 , #1 ) )+BSCALE) ) 

370  XSCREEN=XSCREENADJUST+XSTEP*X 

380  YSCREEN=YSCREENADJUST-Z 

390  IF  TYPE=0  THEN  GOTO  540 

400  IF  Y=0  THEN  GOTO  530 

410  DYSCREEN=YSCREEN-OLDLINE(X, 2) 

420  M=DY  SCREEN/YSTEP 

430  B=YSCREEN-M*XSCREEN 

440  IF  ABS (M)<1  THEN  GOTO  500 

450  IF  DYSCREEN>0  THEN  STEPI=1  ELSE  STEPI=-1 

460  FOR  1=0  TO  CINT(DYSCREEN)  STEP  STEPI 

470  PY=CINT ( OLDLINE (X, 2 )+I) : PX=CINT( (OLDLINE (X, 2 )+I -B)/M) :GOSUB  830 

480  NEXT  I 

490  GOTO  530 

500  FOR  1=1  TO  YSTEP 

510  PX=CINT (OLDLINE (X, 1)+I) : PY=CINT(M* (OLDLINE (X, 1 )+I)+B) :GOSUB  830 

520  NEXT  I 
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530  OLDLINE (X , 1 ) =XSCREEN : OLDLINE (X , 2 ) =YSCREEN 

540  IF  XOFIRSTX  THEN  GOTO  570 

550  PX=CINT(XSCREEN) : PY=CINT (YSCREEN) :G0SUB  830 

560  GOTO  690 

570  DYSCREEN=YSCREEN- YOLDSCREEN 

580  M=DYSCREEN/XSTEP 

590  B=YSCREEN-M*XSCREEN 

600  IF  ABS (M)<1  THEN  GOTO  660 

610  IF  DYSCREEN>0  THEN  STEPI=1  ELSE  STEPI=-1 

620  FOR  1=0  TO  CINT(DYSCREEN)  STEP  STEPI 

630  PY=CINT(YOLDSCREEN+I ) : PX=CINT( (YOLDSCREEN+I-B)/M) :GOSUB  830 

640  NEXT  I 

650  GOTO  690 

660  FOR  1=1  TO  XSTEP 

670  PX=CINT (XOLDSCREEN+I ) : PY=CINT(M*(XOLDSCREEN+I)+B) :GOSUB  830 

680  NEXT  I 

690  XOLDSCREEN=XSCREEN : YOLDSCREEN=YSCREEN 

700  NEXT  X 

710  FOR  ELEMENT=1  TO  639 

720  IF  TEMPBIGGEST ( ELEMENT ) >BIGGEST ( ELEMENT ) 

THEN  BIGGEST ( ELEMENT) =TEMPBIGGEST ( ELEMENT) 

730  IF  TEMPSMALLEST(ELEMENT)<SMALLEST( ELEMENT) 

THEN  SMALLEST ( ELEMENT) =TEMPSMALLEST ( ELEMENT) 

740  TEMPBIGGEST ( ELEMENT )=0 : TEMPSMALLEST( ELEMENT) =1000 

750  NEXT  ELEMENT 

760  NEXT  Y 

770  CLOSE  1 

780  BEEP 

790  A$=INKEY$ : IF  A$=,”,  THEN  790 
800  KILL  " improc#3 . dat" 

810  SYSTEM 
820  END 

830  IF  (PY>BIGGEST(PX) )OR(PY<SMALLEST(PX) ) THEN  PSET(PX,PY) 

840  IF  PY>TEMPBIGGEST ( PX)  THEN  TEMPBIGGEST(PX)=PY 
850  IF  PY<TEMP SMALLEST ( PX)  THEN  TEMPSMALLEST ( PX) =PY 
860  RETURN 
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Appendix  C. 


Hybrid  computer-optical  processing 
with  inexpensive  liquid  crystal  television 

Matt  Young  and  Matt  Weppner 

U.S.  National  Bureau  of  Standards,  Electromagnetic  Technology  Division, 
325  Broadway,  Boulder,  Colorado  80303,  USA 


Abstract 


We  describe  a computer-optical  processing  system  that  uses  an  inexpensive  liquid  crystal 
(LCD)  television  monitor  and  a selective  holographic  filter  for  coherent  pattern  recogni- 
tion. Specifically,  we  use  a digital  computer  to  generate  an  edge  enhanced  image  of  an 
object,  expose  a Fourier  transform  hologram  of  this  image,  and  use  the  hologram  as  a sort  of 
matched  filter  for  recognizing  the  original  object  in  real  time. 

Introduction 


It  may  not  be  true  that  great  minds  think  alike,  but  it  is  demonstrably  true  that  a hand- 
ful of  researchers  have  roughly  simultaneously  discovered  a new  liquid  crystal,  or  LCD, 
television  set  for  use  in  optical  processing. 1-5  The  value  of  such  a device,  of  course,  is 
that  it  is  not  self  luminous  (like  a CRT)  but,  rather,  can  be  used  as  a spatial  light  modu- 
lator for  generating  coherent  images  at  video  rates. 

Typical  LCD  TVs  cost  less  than  US$200  and  may  be  addressed  by  a TV  camera  or  by  a micro- 
computer with  a video  frame  digitizer.  They  therefore  promise  to  allow  coherent  video 
images  and  computer  generated  patterns  for  holography  and  optical  processing  into  almost  any 
optics  laboratory.  We  anticipate  that  they  will  find  their  greatest  application  in  the 
input  plane  of  the  optical  processor,  where  they  are  well  suited,  for  example,  to  matched 
filtering  and  incoherent  correlation  experiments.  Devices  with  higher  resolution  than  the 
currently  available  values  of  about  120  x 140  pixels  (vertical  times  horizontal)  should 
additionally  be  useful  for  change  or  defect  detection6  and,  to  a limited  degree,  to  computer 
generated  holograms.5  In  applications  where  the  transform  plane  can  be  magnified  signifi- 
cantly, the  monitor  may  be  located  there  and  used  as  a computer  generated  spatial  filter  or, 
possibly,  a holographic  filter. 

The  devices  on  the  market  today  are  designed  for  recreation  and  consist  of  a twisted 
nematic  LCD  screen  glued  between  crossed  polarizers.  They  are  viewed  in  transmission  and 
through  a diffuser  inclined  at  an  angle  of  about  45°  to  the  horizon.  Removing  the  diffuser, 
modifying  the  hinge  on  the  device,  and  building  a rigid  frame  takes  several  hours. 

The  LCD  panel  we  use  is  about  55  x 42  mm2  and  uses  a raster  formed  by  a grid  of  fine 
wires;  the  pixels  are  about  0.35  x 0.39  mm2.7  The  polarizers  are  not  flat  and  cause  a 
serious  aberration  in  the  transform  plane1  as  well  as  a lack  of  space  invariance  in  spa- 
tially coherent  systems.4  The  aberration  may  be  corrected  by  contacting  optical  glass 
plates  with  index  matching  fluid  to  both  sides  of  the  display  or  by  removing  the  polarizers 
entirely  and  using  external  polarizers. 

The  fine  wire  raster  causes  many  diffraction  orders  in  the  transform  plane  of  a coherent 
processor;  these  orders  contain  the  information  about  the  raster  itself.  Therefore,  for 
many  applications,  it  will  be  necessary  to  use  low  pass  filtering  to  eliminate  all  but  the 
lowest  diffraction  order  of  the  grid.  For  this  purpose  we  use  a two-stage  system  with  two 
transform  planes.1  The  low  pass  filtering  is  performed  in  the  first  stage;  this  leaves  th._- 
transform  plane  of  the  second  stage  completely  unobstructed  so  that  filters  or  holograph i c 
plates  may  be  easily  located  there. 

The  LCD  display  is  also  well  suited  to  incoherent  correlations,  where  the  object  is  ill 
minated  coherently  but  diffusely,  a hologram  is  recorded,  and  correlations  are  perform'-: 
with  spatially  incoherent  illumination.8'9  (In  thi's  case,  the  quality  of  the  polarizers  is 
irrelevant.)  Such  applications  have  been  hindered  in  the  past  by  the  lack  of  a suital 1<- 
incoherent  display  with  the  same  wavelength  as  the  argon  laser;  the  LCD  monitor  allow:  the 
argon  laser  to  be  used  throughout  the  experiment. 

Hybrid  optical  processor 

The  optical  system  we  use  is  an  extension  of  that  described  in  reference  1 an  1 is 
here  as  Fig.  1.  An  unpolarized  3-mW  He-Ne  laser  beam  is  split  by  an  uncoatei,  w--  1 :«•  : : ■■  . 
splitter;  one  of  the  reflected  beams  is  used  as  the  reference  beam  in  a mat  -he  : f 
experiment.  The  transmitted  beam  is  spatially  filtered  with  a lOx  microscope  ;,«•••♦  i n 
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Figure  1.  Hybrid  processor  consisting  of  a two-stage  Fourier  optical  system  that  includes  a 
liquid  crystal,  or  LCD,  TV  monitor  as  input.  A digital  computer  with  a frame 
digitizer  addresses  the  LCD  monitor  and  a matched  filter  is  recorded  on  the 
photographic  plate.  Camera  TV2  displays  the  cross-correlation  functions  on  a 
monitor  or,  alternatively,  connects  to  the  frame  digitizer,  as  shown  by  the 
dashed  line.  The  three  Fourier  transform  lenses  have  250-mm  focal  length,  f/5.6. 


40-ym  pinhole,  collimated  with  a 360-mm  lens,  and  directed  into  the  two-stage  processor. 

This  system  consists  of  a conventional  4-f  processor  followed  by  a single-lens  processor  in 
which  the  image  is  magnified  slightly  for  ease  of  photographing.  The  lenses  in  the 
processor  are  all  210-mm  copy  or  enlarging  lenses.  A Fresnel  lens  is  inserted  into  the 
final  image  plane  to  serve  as  a field  lens. 

We  located  the  LCD  monitor  in  the  object  plane  of  the  first  stage.  To  optimize  the  sys- 
tem's performance,  we  peeled  the  polarizers  from  the  screen  and  used  external  polarizers. 

We  enclosed  both  the  polarizers  and  the  screen  in  liquid  gates  by  contacting  them  on  both 
sides  to  optical  glass  flats  with  an  index  matching  oil.  The  oil  has  had  no  effect  on  the 
electrical  performance  of  the  monitor. 

The  monitor  is  addressed  by  a microcomputer  with  a video  frame  digitizing  board.  To 
evaluate  the  performance  of  the  monitor,  we  generated  some  patterns,  such  as  that  shown  in 
Fig.  2.  This  pattern  is  radially  symmetric  and  runs  linearly  from  black  in  the  center  to 
white  at  the  edge,  or,  in  terms  of  the  frame  digitizer,  from  0 at  the  center  to  255  at  the 
edge.  (We  used  this  pattern  because  the  monitor  gave  unpredictable  results  with  uniform 
gray  or  with  a gray  scale  whose  gradient  was  either  vertical  or  horizontal.  We  attribute 
this  to  a defect  in  the  video  analog-to-digital  (A-to-D)  converter,  but  we  do  not  know 
whether  or  not  it  is  just  an  idiosyncrasy  of  our  sample.  A newer  model  seems  to  give  some- 
what better  results  in  this  regard.)  After  optimizing  the  polarizers  and  the  monitor's 
brightness  control  and  A-to-D  converters,  we  used  a silicon  detector  to  scan  along  a diam- 
eter of  the  radial  transmittance  pattern.  The  result  is  shown  as  Fig.  3A.  The  monitor  has 
an  overall  contrast  ratio  of  only  10  to  1 and  cannot  respond  to  the  full  dynamic  range  of 
the  TV  camera  or  frame  digitizer.  We  attribute  the  low  contrast  in  part  to  the  presence  of 
the  wires;  these  mainly  add  unwanted  light  to  the  nominally  black  areas,  probably  because  of 
diffraction  or  some  perturbation  of  the  liquid  crystals. 
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Relative  Intensity 


Figure  2.  Radial  transmittance  function  displayed  on  the  LCD  monitor.  A,  before  low  pass 
filtering,  B,  after  low  pass  filtering.  Contrast  ratio  of  B is  5 times  that 
of  A. 


255  0 255 

White  Black  White 


Input  to  LCD  Monitor 


Figure  3.  Transmittance  as  a function  of  position  along  a diameter  of  each  of  tin-  - j 

shown  in  Fig.  2.  Horizontal  scale  is  the  output  of  the  frame  digit  i.-. 
range  from  0 (black)  to  255  (white). 
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When  we  located  a 0.4-mm  low  pass  filter  in  the  transform  plane  of  the  first  stage,  we 
thereby  eliminated  the  wire  grid  from  the  image  (Fig.  2B)  and  measured  a contrast  ratio 
(Fig.  3B)  of  about  50  to  1.  The  monitor,  however,  still  does  not  respond  to  the  full 
dynamic  range  of  the  incoming  TV  signal,  as  is  shown  by  the  flat  portions  of  the  curve  near 
the  center  and  the  edges  of  the  display.  These  signify  a relative  lack  of  both  shadow  and 
highlight  detail;  adjusting  the  brightness  or  the  video  A-to-D  converter  level  simply  shifts 
the  sloped  portions  of  the  curve  horizontally  toward  or  away  from  the  center  of  the  graph. 

To  complete  the  system,  we  located  a holographic  plate  in  the  transform  plane  of  the 
second  stage  and  introduced  the  reference  beam  through  a pair  of  polarizers  for  intensity 
control.  (The  reference  beam  angle  of  about  15°  was  determined  by  the  size  of  the  transform 
lens.)  For  these  experiments,  we  did  not  collimate  the  reference  beam,  so  it  had  a waist  at 
the  location  of  the  laser,  about  2 m from  the  hologram.  This  may  reduce  the  position 
invariance  of  the  system,  but  we  did  not  check  for  position  invariance. 

Experiment 

For  our  experiment,  we  illuminated  a pair  of  pliers  diffusely  with  two  incandescent  lamps 
and  white,  matte  paper  diffusers.  Using  a TV  camera  fitted  with  a zoom  lens,  we  captured 
and  stored  an  image  of  the  pliers  with  the  computer  and  the  frame  digitizer.  The  image  dis- 
played on  the  LCD  monitor  is  shown  in  Fig.  4A  prior  to  low  pass  filtering  with  the  optical 
system.  Figure  4B  shows  the  same  image  after  spatially  filtering  with  the  0.4-mm  aperture 
to  eliminate  the  wire  grid.  When  the  low  pass  filter  is  chosen  properly,  only  the  grid  is 
eliminated  from  the  picture;  there  is  no  loss  of  resolution  of  the  image  itself. 


Figure  4.  Pair  of  pliers  used  as  a subject  for  pattern  recognition  experiment  and  displayed 
on  the  LCD  monitor.  A,  image  of  the  pliers  themselves.  B,  low  pass  filtered  to 
remove  the  pixels  on  the  LCD  screen.  C,  binary  image  of  the  pliers,  wherein 
values  below  a threshold  are  set  equal  to  0 and  values  above  that  threshold  are 
set  equal  to  255.  D,  edge  enhanced  model  of  the  pliers. 
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Next  we  used  the  computer  to  prepare  a clipped,  or  binary,  image  of  the  pliers,  Fig.  4C, 
where  intensities  above  a certain  level  are  set  equal  to  255  (white)  and  those  below  that 
level  are  set  equal  to  0 (black).  In  the  one  dimensional  analogy, 


Next,  we  performed  a two  dimensional  edge  enhancement  in  which  we  retained  only  the 
inside  edge  of  the  binary  image;  that  is,  the  image  is  black  except  for  a white  region  that 
runs  along  the  inside  edge  of  the  binary  image.  Symbolically, 


(2) 


where  the  asterisk  (*)  denotes  cross  correlation  and  values  less  than  zero  are  interpreted 
as  zero.  We  refer  to  the  resulting  image  as  the  model;  it  is  shown  in  Fig.  4D.  The  width 

of  the  enhanced  edge  is  determined  by  the  number  of  zeros  between  the  peaks  of  the  kernel  of 

the  cross  correlation  integral.  We  chose  this  width  to  be  two  pixels  on  the  LCD  monitor; 
this  corresponds  to  about  six  pixels  on  a conventional  monitor  with  the  full  resolution  of 
about  380  pixels. 

We  used  the  model,  somewhat  in  the  manner  of  references  8 and  9,  as  the  object  in  a pat- 
tern recognition  experiment.  We  inserted  a holographic  plate  in  the  second  transform  plane 

and  recorded  a Fourier  transform  hologram  of  the  model.  The  value  of  using  the  model, 

rather  than  the  object  itself,  is  that  the  cross  correlation  function  of  the  model  and  the 
object , 


(3) 


is,  with  an  irregular  two  dimensional  object,  sharper  than  the  autocorrelation  function  of 
the  object  itself. 


(4) 


and  permits  greater  discrimination  against  similar  objects.  In  addition,  the  model  is 
located  on  a black  background,  and  its  Fourier  transform  has  a relatively  weak  dc  or  zer 
order  component.  As  a result,  the  exposure  of  the  photographic  plate  is  comparatively  mi- 
form  and  allows  a hologram  to  be  recorded  with  roughly  constant  diffraction  efficiency 
all  relevant  spatial  frequencies. 


Results 


Placing  the  emulsion  side  away  from  the  transform  lens,  we  exposed  some  Fourier  • : 
holograms  to  serve  as  matched  filters  for  the  model.  We  processed  some  of  the-  it.  • ' u- 

sium  ferricyanide  bleach  (15  mg/mL  for  4 min).  These  may  be  expected  to  have  > !i  f ! : .>  • 
efficiency  perhaps  10  times  that  of  amplitude  transmission  holograms  and  ar<-  :•  i 
tolerant  of  variations  of  exposure.10  Our  best  filter  was  just  such  a phase  h u i 
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Using  a second  TV  camera  denoted  TV2  and  shown  in  Fig.  1,  we  viewed  the  autocorrelation 
function  of  the  model  on  a monitor.  Because  the  image  is  mostly  black,  we  switched  the 
camera's  automatic  gain  control  off  to  avoid  overexposing  the  autocorrelation  spot.  The  AGC 
was  left  off  during  all  the  observations  that  follow.  We  did  not  need  to  synchronize  the 
camera  TV2  to  the  computer,  evidently  because  the  decay  time  of  the  LCD  panel  was  longer 
than  the  frame  rate  of  1/60  s. 

Figure  5A  shows  a photograph  of  the  central  portion  of  the  monitor.  The  spot  is  sharp 
and  has  a width  of  about  six  TV  lines;  this  agrees  well  with  the  expected  value  of  two 
pixels  on  the  lower  resolution  LCD  monitor.  The  spot  can  also  be  seen  clearly  with  the  eye 
and  is  surrounded  by  a diffraction  ring  that  does  not  appear  in  the  photograph. 

Next,  we  displayed  the  clipped  image  of  the  object  on  the  LCD  monitor,  without  moving  the 
apparatus  or  adjusting  either  TV2  or  the  brightness  of  the  monitor.  The  resulting  cross 
correlation  function  is  shown  in  Fig.  5B.  We  then  switched  to  real  time  and  displayed  the 
real  object  on  the  LCD  monitor.  The  relative  brightness  of  the  object  was  substantially 
less  than  the  255,  or  white,  of  the  computer  processed  images.  The  cross  correlation  with 


Figure  5.  Correlation  functions  as  seen  by  TV2  displayed  on  a monitor.  A,  autocorrelation 
function  of  the  model.  B-D,  cross  correlation  function  of  the  model  with  (B) 
clipped  tool,  (C)  tool  in  real  time,  and  (D)  tool  rotated  2°.  All  photographs 
were  taken  with  the  same  exposure. 


Figure  6.  Intensity  as  a function  of  position  along  a horizontal  TV  line  through  the  center 
of  the  cross  correlation  peak  of  Fig.  5C.  Average  of  four  frames. 
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the  model  was  weaker  but  plainly  visible  on  the  monitor  (Fig.  5C).  Finally,  we  rotated  the 
object  2°.  The  cross  correlation  function  in  this  case  was  still  visible  and  is  shown  in 
Fig.  5D. 

We  then  switched  the  camera,  TV2,  to  the  computer  and  digitized  its  output.  We  could  not 
record  the  autocorrelation  function  of  the  model,  because  that  would  have  required  two  frame 
digitizers;  indeed,  this  is  one  of  the  reasons  that  we  designed  our  experiment  to  examine 
the  cross  correlation  of  the  unprocessed  object,  rather  than  a processed  image,  with  the 
model . 

We  have  made  no  attempt  to  optimize  the  diffraction  efficiency  of  the  hologram.  What  was 
plainly  visible  on  the  monitor  was  not  so  plain  to  the  computer;  the  eye  does  a marvelous 
job  of  smoothing  and  averaging.  We  therefore  digitized  four  frames  and  averaged  them.  Fig- 
ure 6 is  a graph  of  the  relative  intensity  across  one  horizontal  line  in  the  averaged  pic- 
ture, multiplied  by  3 for  clarity  of  display.  The  constant  value  of  about  115  is  simply 
dark  current  from  the  TV  camera. 

Figure  7 shows  a histogram  of  the  entire  screen  (of  which  Fig.  6 is  a single  line).  (All 
values  greater  than  zero  are  shown  as  small  ticks,  whether  or  not  they  properly  round  to 
zero.)  The  highest  value,  228,  represents  the  cross  correlation  peak.  The  average  value  of 
the  dark  current  is  about  118,  and  the  dark  noise  fluctuations  range  from  roughly  104  to 
132.  If  we  take  this  to  be  ±3  sigma,  we  estimate  the  signal-to-noise  ratio  (SNR)  to  be 
approximately  24. 


Figure  7.  Histogram  of  the  entire  screen  from  which  Fig.  6 is  derived.  The  signal  peat  is 
228,  the  noise  average  118,  and  ±3  sigma  about  28.  The  signal-to-noise  rati  is 
about  24. 


Discussion 


We  have  used  a digital  computer  to  generate  an  edge  enhanced  (but  non-negative)  !el  i 
a tool,  prepared  a Fourier  transform  hologram  of  the  model,  and  cross  correlated  the  ie 
in  real  time  with  an  image  of  the  original  tool. 

We  might  have  achieved  better  discrimination  if  we  had  either  edge  enhance  1 • he  • 
real  time  (as  in  reference  8)  or  edge  enhanced  with  positive  and  negative  val  its  i- 
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reference  9).  For  practical  reasons,  we  did  neither.  The  latter,  in  particular,  would 
have  required  an  encoding  scheme  that  required  a space-bandwidth  product  greater  than  the 
LCD  monitor  offered  and  would  not  have  achieved  significantly  different  results  unless  we 
had  used  the  equivalent  of  heterodyne  detection.  We  could  have  achieved  a similar  result  by 
edge  enhancing  optically  (with  a bandpass  filter  in  the  first  transform  plane);  this 
achieves  a phase  shift  in  the  electric  field  at  the  location  of  an  edge  (as  is  shown  by  the 
presence  of  a zero  of  intensity  at  the  edge11).  We  chose,  however,  to  edge  enhance  by  com- 
puter since  we  thought  that  this  would  be  more  readily  controllable  and  result  in  a simpler, 
singly  peaked  intensity  distribution  at  the  edge. 

Finally,  we  have  performed  our  experiment  using  coherent  processing.  This  is  far  more 
efficient  than  incoherent  processing,  which  requires  diffused  illumination,  and  has  allowed 
us  to  use  a small  He-Ne  laser  rather  -than  the  much  higher  power  argon  laser  used  in  refer- 
ences 8 and  9.  Except  for  the  limited  space-bandwidth  product  offered  by  the  LCD  monitor, 
most  other  considerations  are  about  equal,  and  the  use  of  the  small  and  inexpensive  He-Ne 
laser  may  therefore  militate  in  favor  of  coherent  processing  for  this  sort  of  pattern  recog- 
nition . 
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